#  Building models of quarks and gluons with an arbitrary number of colors
using Cartan-Polyakov loops

Renan Câmara Pereira<EMAIL_ADDRESS>Centro de Física da
Universidade de Coimbra (CFisUC), Department of Physics, University of
Coimbra, P-3004 - 516 Coimbra, Portugal  Pedro Costa<EMAIL_ADDRESS>Centro de
Física da Universidade de Coimbra (CFisUC), Department of Physics, University
of Coimbra, P-3004 - 516 Coimbra, Portugal

###### Abstract

In this work we introduce the concept of Cartan-Polyakov loops, a special
subset of Polyakov loops in the fundamental and antifundamental representation
of the $\mathrm{SU}(N_{c})$ group, $\ell_{\mathrm{F},k}$ and
$\overline{\ell}_{\mathrm{F},k}$ respectively, with charges
$k=1,\ldots,(N_{c}-1)/2$. It constitutes a sufficient set of independent
degrees of freedom and it is used to parametrize the thermal Wilson line.
Polyakov loops not contained in this set are classified as non-Cartan-Polyakov
loops. Using properties of the characteristic polynomial of the thermal Wilson
line, we write a non-Cartan-Polyakov loop charge decomposition formula. This
formalism allows one to readily build effective models of quarks and gluons
with an arbitrary number of colors. We apply it to the
Polyakov$-$Nambu$-$Jona-Lasinio model and to an effective glue model, in the
mean field approximation, showing how to directly extend these models to
higher values of $N_{c}$.

##  I Introduction

Gauge theories with the symmetry group $\mathrm{SU}(N)$, are fundamental
building blocks to our current understanding of the laws of physics at the
microscopic level. The theory of strong interactions, Quantum Chromodynamics
(QCD), is a perfect example of a gauge theory based on the local invariance
under transformations of the $\mathrm{SU}(N_{c})$ symmetry group, with $N_{c}$
being the number of colors. This theory has several important properties like
scaling, asymptotic freedom, chiral symmetry breaking and confinement [1].

The interplay of the properties of QCD alongside its non-perturbative nature
at low energies, leads the way to an extremely rich phase structure for
strongly interacting matter, also known as the QCD phase diagram [2].
Originally, only confined and deconfined phases were displayed in the diagram
however, since then, other features have been discovered or proposed, such as,
the existence of the quark-gluon plasma at high temperatures, the existence of
the conjectured critical endpoint at finite temperature and baryon density,
the possible diquark condensation in high density baryon matter or the
confined chirally symmetric phase [3], to name a few.

However, in order to study QCD in these energy regimes we are limited in
options due to the non-perturbative behavior of the theory. From a theoretical
point of view, one can use lattice QCD calculations, Dyson$-$Schwinger
equations and effective models. Lattice methods are the main tool to
understand strongly interacting matter from first principles, at least for
finite temperature and zero chemical potential. At finite density, lattice QCD
suffers from the infamous sign problem which renders the importance sampling,
necessary in Monte Carlo simulations, ineffective [4].

Deviating from the real world of three colors and studying systems with
increasing $N_{c}$, or considering the large $N_{c}$ limit (also known as ’t
Hooft limit, originally introduced in Ref. [5]), has been demonstrated to be a
really good tool to better understand the behavior of the QCD properties [6,
7, 8, 9, 10, 11, 12], allowing for the exploration of the non-perturbative
behavior of QCD in a systematic way [6, 7, 9]. The large $N_{c}$ limit is
taken by considering a large number of colors, while fixing
$\lambda=g^{2}N_{c}$, with $g$ the strong coupling and $\lambda$ the so-called
’t Hooft coupling [5].

Indeed, the large $N_{c}$ limit of the theory should still exhibit confinement
and chiral symmetry breaking since the theory remains asymptotically free,
with a growing coupling at low energies [7, 13]. Not only that, some lattice
studies, the phenomenological OZI (Okubo$-$Zweig$-$Iizuka) rule and the
reproduction of ratios of low-lying meson masses [14], also seem to point in
the direction of a large $N_{c}$ limit of QCD resembling the real world case
of $N_{c}=3$ [15, 16, 17, 18, 19, 13]. Although QCD in the large $N_{c}$ limit
is a simpler theory, even in such approximation it has not yet to been solved
[20, 21]. Non-abelian gauge theories, with the gauge group $\mathrm{SU}(N)$,
in the large $N$ limit are also very important in the study of string
theories, in particular, in studies involving the conjectured AdS/CFT
correspondence, i.e., the duality between gauge theories and theories of
gravity [22, 23, 9]. Indeed, one can also use this correspondence to study the
QCD phase diagram [24].

In the spirit of the large $N_{c}$ limit, Pisarski et al. proposed the
existence of a new phase of matter in the phase diagram of strongly
interacting matter, the so-called quarkyonic phase, expected when the number
of colors, $N_{c}$, is very large [25]. It was conjectured that, in this
limit, the QCD phase diagram can be separated in three distinct phases, two at
lower temperatures and one at high temperatures. The latter is the usual
quark-gluon plasma, constituting the deconfined phase of matter. At low
temperatures there would be two distinct phases: a mesonic phase, described by
confined matter with zero baryon number density and, at higher chemical
potentials, the so-called quarkyonic phase, described by confined matter with
finite baryon number density [25, 26]. The conjectured quarkyonic phase of
matter was discussed and studied in Refs. [26, 27, 28, 24, 29], being also
explored in the context of neutron stars in Refs. [30, 31].

From the phenomenological point of view, several effective models of QCD have
been built in order to study the properties of the conjectured quarkyonic
phase. Some of these models were based on the Polyakov$-$Nambu$-$Jona-Lasinio
(PNJL) model, in which quark-quark interactions are point-like, with a static
gluon field in the temporal direction acting as a background field. These
models possess all the global symmetries of QCD, including the
confinement/deconfinement transition in the sense of the breaking of the
$\mathrm{Z}(N_{c})$ symmetry and which can be measured by the expectation
value of the Polyakov loop. This quantity is an exact order parameter in the
pure glue theory but only approximate when quarks are introduced in the
Lagrangian. The Polyakov loop is extremely used as an order parameter to study
the deconfinement transition in both ab-initio calculations like lattice QCD,
in Dyson-Schwinger equations, and in phenomenological models like the
Polyakov-quark-meson model or the PNJL model.

In this work we show how to build models of quarks and gluons with an
arbitrary number of colors, in a systematic way. More precisely, we will build
a version of the Polyakov$-$Nambu$-$Jona-Lasinio model, in the mean field
approximation and consider a glue effective potential, both for any number of
colors. Being able to study such models systematically can be quite useful in
order to better understand the deconfinement phase transition, the conjectured
quarkyonic phase and the properties of quarks and gluons.

##  II Center symmetry and the Polyakov loop

In the imaginary time formalism, often used to study quantum field theory at
finite temperature, boson and fermion fields must respect periodic and
antiperiodic boundary conditions, respectively. Using this prescription, it
can be shown that, at the Lagrangian level, the boson fields of the Yang-Mills
theory have symmetric boundary conditions with respect to transformations of
the discrete group $\mathrm{Z}(N_{c})$, the center of the
$\mathrm{SU}(N_{c})$. The introduction of fermion fields, on the other hand,
explicitly breaks such invariance. Hence, Yang-Mills theory at finite
temperature have boundary conditions that are invariant under
$\mathrm{Z}(N_{c})$ transformations while QCD does not share the same
property, due to the presence of the quark fields. However, even in the pure
glue theory, such symmetry can be spontaneously broken.

The Polyakov loop, $\Phi$, is a gauge invariant quantity, that can be used as
an exact order parameter for the spontaneous symmetry breaking of the
$\mathrm{Z}(N_{c})$ symmetry in the Yang-Mills theory [32]: if the Polyakov
loop is zero, the $\mathrm{Z}(N_{c})$ symmetry is preserved while, if it is
non-zero, the symmetry is broken.

One can check the gauge invariance of the Polyakov loop by considering the
gauge transformation of the thermal Wilson loop, $L\quantity(\bm{x})$, from
which one can define the Polyakov loop. The Wilson line can be defined along a
path in spacetime, from an initial position ($x_{i}$) to a final position
($x_{f}$), as follows:

$\displaystyle
W(x_{i},x_{f})=\mathcal{P}\exp\quantity[ig\int_{x_{i}}^{x_{f}}\differential{z^{\mu}}A_{\mu}(z)].$
(1)

Here, $\mathcal{P}$ is the path-ordering operator, $g$ is the strong
interaction coupling and $A_{\mu}$ is the gauge field111 Under a gauge
transformation, the gauge field $A_{\mu}(x)$ transforms as:  $\displaystyle
A_{\mu}(x)\xrightarrow[]{\text{g.t.}}A_{\mu}^{\prime}(x)=\Omega(x)\quantity(A_{\mu}+\frac{1}{ig}\partial_{\mu})\Omega^{\dagger}(x).$
(2) . The Wilson line of the gauge transformed field is given by:

$\displaystyle
W(x_{i},x_{f})\xrightarrow[]{\text{g.t.}}W^{\prime}(x_{i},x_{f})=\Omega(x_{f})W(x_{i},x_{f})\Omega^{\dagger}(x_{i}).$

By considering a closed path ($x_{i}=x_{f}=x$), one gets the so-called Wilson
loop, which is a gauge invariant quantity due to the cyclic property of the
trace:

$\displaystyle\tr\quantity[W(x_{i},x_{f})]\xrightarrow[]{\text{g.t.}}\tr\quantity[W^{\prime}(x_{i},x_{f})]$
$\displaystyle=\tr\quantity[\Omega(x_{f})W(x_{i},x_{f})\Omega^{\dagger}(x_{i})]$
$\displaystyle=\tr\quantity[W(x_{i},x_{f})].$ (3)

Considering the Wilson loop wrapped in imaginary time (i.e., at finite
temperature with periodic boundary conditions) one defines the thermal Wilson
line,

$\displaystyle
L\quantity(\bm{x})=\mathcal{P}\exp\quantity[i\int_{0}^{\beta}\differential{\tau}A_{4}(\bm{x},\tau)].$
(4)

Here, $A_{4}=ig{A}_{\mu}^{a}T_{a}\delta_{0}^{\mu}$, is the gluon field in the
temporal direction ($T_{a}$ are the $\mathrm{SU}(N_{c})$ generators in the
fundamental representation, ${A}_{\mu}^{a}$ are the gluon fields in the
adjoint representation and the index $a$ runs from $0$ to $N_{c}^{2}-1$). This
quantity is also gauge invariant. Finally, one can define the Polyakov loop on
a spatial coordinate $\bm{x}$, in the fundamental representation of the
$\mathrm{SU}(N_{c})$ group, $\Phi(\bm{x})$, as the trace of the thermal Wilson
loop with the field $A_{4}$ written in the fundamental representation [32,
33]:

$\displaystyle\Phi(\bm{x})=\frac{1}{N_{c}}\tr\quantity[L_{\mathrm{F}}\quantity(\bm{x})].$
(5)

Here, the trace is made over color space and the index $\mathrm{F}$ indicates
that the matrix $A_{4}$ is written in the fundamental representation of the
$\mathrm{SU}(N_{c})$ group. The Polyakov loop transforms akin to the fermion
fields and only the identity element of the $\mathrm{Z}(N_{c})$ group
maintains the symmetry indeed, under a global $\mathrm{Z}(N_{c})$
transformation, it transforms as a field with charge one [32, 34].

The physical interpretation behind the two previously discussed phases,
corresponding to zero and non-zero Polyakov loop, can be understood by
evaluating the Helmholtz free energy of a single quark in a background of
gluons. To this end, one can solve the static Dirac equation coupled to a
static background field and use it to calculate the Helmholtz free energy of a
single quark in a background of gluons [35, 32, 36]. It yields that, if the
Polyakov loop is finite, a single quark can be created in the gluonic
background with a finite amount of energy while, if the Polyakov loop is zero,
an infinite amount of energy is required [32, 36]. Hence, a finite Polyakov
loop would correspond to a deconfined phase of matter while, a zero Polyakov
loop to a confined phase.

As already pointed out, the Polyakov loop is often incorporated in
phenomenological models of QCD, like the PNJL or quark-meson models, in order
to study the properties of the phase diagram of strongly interacting matter
and explore the relation between the restoration of chiral symmetry and the
deconfinement transition. In the following sections we will show how to build
both a PNJL model and an effective glue potential, with an arbitrary number of
colors, which can be used to explore the properties of strongly interacting
matter.

##  III Cartan-Polyakov loops

In this section, we will introduce the concept of Cartan-Polyakov loops which
will enable one to readily derive quark and gluon models with an arbitrary
number of colors.

Consider the thermal Wilson line defined in Eq. (4), $L$. Allying the unity
determinant property of this quantity (it is an element of the
$\mathrm{SU}(N_{c})$ group) with the fact that, at every spacetime point, this
matrix can be diagonalized with a gauge transformation [32], one can conclude
that this quantity has a total of $N_{c}-1$ degrees of freedom, the same
number as the rank of the $\mathrm{SU}(N_{c})$ group [37]. Knowing the value
of $N_{c}-1$ traces of different powers of this matrix, allows one to build a
system of $N_{c}-1$ equations to the same number of unknown matrix entries.
Hence, the traces can be used to parametrize the thermal Wilson line [38].

As introduced earlier, the Wilson line in the fundamental representation
transforms as a field with charge $+1$, under a $\mathrm{Z}(N_{c})$ symmetry
transformation. The antifundamental representation transforms as a field with
charge $-1$. As pointed out in Ref. [39], the charges are defined modulo
$N_{c}$, meaning that the antifundamental charge, $-1$, is equivalent to the
charge $N_{c}-1$.

One can build Polyakov loops in higher dimensional representations, using
combinations of fundamental and antifundamental Wilson lines. Examples include
Wilson lines in the adjoint representation of $\mathrm{SU}(N_{c})$, two-index
symmetric, antisymmetric and other representations, with higher index counts.
These different possible representations transform under $\mathrm{Z}(N_{c})$
as fields with different charges. For more information on different
representations of the Wilson line, see [39, 40, 41]. In particular, in Ref.
[39], it was shown how the classification procedure can be made in the large
$N_{c}$ limit.

In $\mathrm{SU}(3)$, the fundamental and antifundamental Polyakov loops
constitute a sufficient set from which one can build any other higher
representation Polyakov loop. For larger values of $N_{c}$ this is not the
case [42] and a higher number of degrees of freedom must be considered. In
this work we will use Polyakov loops with different charges, in the
fundamental representation, in order to easily build effective models
containing both gluon and quark degrees of freedom, for an arbitrary $N_{c}$.
More specifically, we will show how to build the Polyakov$-$Nambu$-$Jona-
Lasinio model, in the mean field approximation and an effective gluon
potential, for increasing values of $N_{c}$. Such constructions can be very
helpful in better understanding the thermodynamics of QCD at large $N_{c}$, as
well as, studying the properties of the conjectured quarkyonic phase in the
QCD phase diagram or determining the behavior of the critical endpoint for
increasing $N_{c}$ [29].

We define the Polyakov loop with charge $k$, as the trace of the $k-$power of
the thermal Wilson line, in a particular representation, $\mathcal{R}$, of the
$\mathrm{SU}(N_{c})$ group. Polyakov loops with negative charge, are defined
analogously, as the trace of the $k-$power of the adjoint thermal Wilson line
[34]. Formally, we write:

$\displaystyle\ell_{\mathcal{R},k}\quantity(\bm{x})=\frac{1}{N_{\mathcal{R}}}\tr[L_{\mathcal{R}}(\bm{x})^{k}\big{]},$
(6)
$\displaystyle\overline{\ell}_{\mathcal{R},k}\quantity(\bm{x})=\frac{1}{N_{\mathcal{R}}}\tr[L_{\mathcal{R}}^{\dagger}(\bm{x})^{k}\big{]}.$
(7)

Here, $N_{\mathcal{R}}$ is the dimension of the representation. In the
literature, the Polyakov loop variable in the fundamental representation is
denoted by $\Phi$ and $\overline{\Phi}$, see Eq.(5). Thus, in our notation,
the following equivalence holds: $\Phi=\ell_{\mathrm{F},1}$ and
$\overline{\Phi}=\overline{\ell}_{\mathrm{F},1}$. The zero charge Polyakov
loop is trivially defined by the normalized trace of the
$N_{\mathcal{R}}\times N_{\mathcal{R}}$ identity matrix,
$\ell_{\mathcal{R},0}=\overline{\ell}_{\mathcal{R},0}=1$. In other words, the
Polyakov loop with charge $k$ is the normalized trace of a thermal Wilson line
wrapped $k$ times around the compact imaginary time dimension [20]:

$\displaystyle\ell_{\mathcal{R},k}\quantity(\bm{x})=\frac{1}{N_{\mathcal{R}}}\tr\mathcal{P}\exp\quantity[i\int_{0}^{k\beta}\differential{\tau}A_{4}(\bm{x},\tau)].$
(8)

Interestingly, for a given number of colors, these objects are not independent
of each other. Indeed, a given Polyakov loop with charge $k$ can be written as
a function of other Polyakov loops, with different charges. This is connected
to a previously mentioned property of the thermal Wilson line: it has a total
of $N_{c}-1$ independent degrees of freedom. As discussed earlier, one can
parametrize the thermal Wilson line by calculating the trace of $N_{c}-1$
powers of this matrix. Since a Polyakov loop with charge $k$ is simply the
trace of the $k-$power of the Wilson thermal line, there is a total of
$N_{c}-1$ independent charges. This immediately poses the question about which
charges should be chosen as degrees of freedom to parametrize a
phenomenological model.

This motivates the introduction of a special subset of Polyakov loops, which
we denote by Cartan-Polyakov loops. For a specific number of colors, $N_{c}$,
we denote the subset of $N_{c}-1$ Polyakov loops, with charges chosen to be
independent of each other, by Cartan-Polyakov loops and the remaining ones, by
non-Cartan-Polyakov loops. In this work we define the Cartan-Polyakov loop set
to be constituted by Polyakov loops in the fundamental and antifundamental
representations of the $\mathrm{SU}(N_{c})$ group, $\ell_{\mathrm{F},k}$ and
$\overline{\ell}_{\mathrm{F},k}$, with charges $k=1,\ldots,(N_{c}-1)/2$.
Hence, the Cartan-Polyakov loop set, $\mathcal{C}_{N_{c}}$, is given by:

$\displaystyle\mathcal{C}_{N_{c}}=\quantity{\ell_{\mathrm{F},1},\ldots,\ell_{\mathrm{F},(N_{c}-1)/2},\overline{\ell}_{\mathrm{F},1},\ldots,\overline{\ell}_{\mathrm{F},(N_{c}-1)/2}}.$
(9)

Although this choice is made a priori and may seem arbitrary, it is guided by
simplicity: although other choices exist222For instance, one could choose
$\ell_{\mathrm{F},k}$ and $\overline{\ell}_{\mathrm{F},k}$ with
$k=(N_{c}-1)/2,\ldots,(N_{c}-1)$ or, only consider the fundamental
representation, $\ell_{\mathrm{F},k}$ with $k=1,\ldots,(N_{c}-1)$., this
particular one produces straightforward one-to-one mappings from the Cartan
loops to non-Cartan loops while, for other choices, this may not be the case.
For $N_{c}=3$, for example, the Cartan loops, are given by the set
$\mathcal{C}_{3}=\quantity{\ell_{\mathrm{F},1},\overline{\ell}_{\mathrm{F},1}}$.
For the cases in which $N_{c}$ is an even number, we opted to remove the
element $\overline{\ell}_{\mathrm{F},(N_{c}-1)/2}$ from the set defined in Eq.
(9). The reason for this is based on the fact that, for an even number of
colors, there is freedom to choose between $\ell_{\mathrm{F},(N_{c}-1)/2}$ and
$\overline{\ell}_{\mathrm{F},(N_{c}-1)/2}$. In other words, picking one of
them to include in the Cartan set allows one to write the other in terms of
the other loops in the Cartan set. In this work, for simplicity, we mainly
focus on results for odd $N_{c}$: for such case there is a symmetry between
the number of $\ell_{\mathrm{F},k}$ and $\overline{\ell}_{\mathrm{F},k}$
contained in the Cartan set.

In the following, we will write recurrence relations, which will allow us to,
within a fixed representation of the $\mathrm{SU}(N_{c})$ group, decompose a
non-Cartan-Polyakov loop with charge $k$, in terms of Polyakov loops with
smaller charges, $k-1,k-2,\ldots$ The idea is to use such recurrence relation
in order to build effective models of quarks and gluons whose only degrees of
freedom are the Cartan-Polyakov loops. For example, when building a model with
$N_{c}=5$, the Cartan loops, are given by the set
$\mathcal{C}_{5}=\quantity{\ell_{\mathrm{F},1},\ell_{\mathrm{F},2},\overline{\ell}_{\mathrm{F},1},\overline{\ell}_{\mathrm{F},2}}$.
However, in the construction process, non-Cartan-Polyakov loops such as,
$\ell_{\mathrm{F},3}$, $\ell_{\mathrm{F},4}$,
$\overline{\ell}_{\mathrm{F},3}$, $\overline{\ell}_{\mathrm{F},4}$, will
arise. Using the recurrence relation that we will derive in the next section,
enables one to write these non-Cartan-Polyakov loops as functions of loops
with smaller charge, more specifically the ones present in the set
$\mathcal{C}_{5}$.

We note that it may not be always possible to write effective models in terms
of Cartan-Polyakov loops. Indeed in some models arises the necessity to
involve elements of the thermal Wilson line matrix [41]. The mapping between
the matrix elements of $L_{\mathcal{R}}$ and the Cartan-Polyakov loops is not
straightforward.

###  III.1 Decomposition of non-Cartan-Polyakov loops

In order to derive the non-Cartan-Polyakov loop charge decomposition formula,
we will make use of an identity between the coefficients of the characteristic
polynomial of an unitary matrix. Consider the characteristic polynomial of the
thermal Wilson line in a particular representation, $L_{\mathcal{R}}$. It can
be written as (see Appendix A):

$\displaystyle p\quantity(L_{\mathcal{R}},\lambda)$
$\displaystyle=\det\quantity[L_{\mathcal{R}}-\lambda I]$
$\displaystyle=\quantity(-1)^{N_{\mathcal{R}}}\quantity(\lambda^{N_{\mathcal{R}}}+\sum_{i=1}^{{N_{\mathcal{R}}}-1}c_{i}(L_{\mathcal{R}})\lambda^{{N_{\mathcal{R}}}-i}+\quantity(-1)^{N_{\mathcal{R}}}\det\quantity[L_{\mathcal{R}}]).$
(10)

Here, $c_{i}(L_{\mathcal{R}})$ are coefficients calculated using (see Eq.
(LABEL:c_k_definition)):

$\displaystyle c_{k}(L_{\mathcal{R}})$
$\displaystyle=\frac{(-1)^{k}}{k!}\det\quantity[J_{\mathcal{R},k}].$ (11)

The matrix $J_{\mathcal{R},k}$ will be essential in the derivation of the
reducing formula. It is a $k\times k$, square matrix defined by (see Eq.
(LABEL:c_k_definition)):

$\ell_{\mathcal{R},k}$ was defined in Eq. (6). One can also write it in terms
of matrix elements, $(J_{\mathcal{R},k})_{i,j}$ as,

$\displaystyle(J_{\mathcal{R},k})_{i,j}=N_{\mathcal{R}}\ell_{\mathcal{R},i-j+1}H\quantity(i-j)+i\delta_{i,j-1}.$
(13)

Here, $H\quantity(n)$ is the discrete version of the Heaviside step function,
$H\quantity(n)=0$ if $n<0$ or $H\quantity(n)=1$ if $n\geq 0$.

We will also need the matrix $\overline{J}_{\mathcal{R},k}$, which has an
identical structure to the one defined above but with the replacement,
$\ell_{\mathcal{R},k}\to\overline{\ell}_{\mathcal{R},k}$. Its definition comes
from the characteristic polynomial of the adjoint thermal Wilson line,
$p(L_{\mathcal{R}}^{\dagger},\lambda)$, with coefficients given by,
$c_{k}(L_{\mathcal{R}}^{\dagger})=\nicefrac{{(-1)^{k}}}{{k!}}\det\quantity[\overline{J}_{\mathcal{R},k}]$.

Excluding the upper triangular section of the matrix $J_{\mathcal{R},k}$, all
its entries are Polyakov loops in the representation $\mathcal{R}$, with
different charges, ranging from $1$ to $k$. Also, the higher charge Polyakov
loop, $\ell_{\mathcal{R},k}$, is only present in the entry $(k,1)$. All the
other entries in the matrix are numbers or Polyakov loops of charge smaller
than $k$. The idea is to isolate this quantity from the matrix
$J_{\mathcal{R},k}$ and derive a closed expression for $\ell_{\mathcal{R},k}$
in terms of Polyakov loops of smaller charge i.e.,
$\ell_{\mathcal{R},k-1},\ell_{\mathcal{R},k-2},\ldots,\ell_{\mathcal{R},1}$
and
$\overline{\ell}_{\mathcal{R},k-1},\overline{\ell}_{\mathcal{R},k-2},\ldots,\overline{\ell}_{\mathcal{R},1}$.

To this end, consider the Laplace expansion of the determinant on the first
column of the matrix $J_{\mathcal{R},k}$. It can be written as:

$\displaystyle\det\quantity[J_{\mathcal{R},k}]$
$\displaystyle=\sum_{i=1}^{k}\quantity(-1)^{i+1}N_{\mathcal{R}}\ell_{\mathcal{R},i}\det\quantity[J_{\mathcal{R},k}^{(i,1)}]$
$\displaystyle=\quantity(-1)^{k+1}N_{\mathcal{R}}\ell_{\mathcal{R},k}\det\quantity[J_{\mathcal{R},k}^{(k,1)}]+N_{\mathcal{R}}\sum_{i=1}^{k-1}\quantity(-1)^{i+1}\ell_{\mathcal{R},i}\det\quantity[J_{\mathcal{R},k}^{(i,1)}].$
(14)

Here, $J_{\mathcal{R},k}^{(i,j)}$ is the square $(k-1)\times(k-1)$ matrix,
obtained by removing the $i-$th line and $j-$th column of matrix
$J_{\mathcal{R},k}$. The determinant in the first term,
$\det[J_{\mathcal{R},k}^{(k,1)}]$, can be readily calculated by the product of
the diagonal entries since the matrix $J_{\mathcal{R},k}^{(k,1)}$ is of the
lower triangular type. It yields:

$\displaystyle\det\quantity[J_{\mathcal{R},k}^{(k,1)}]=1\times 2\times
3\times\ldots\times(k-1)=\frac{k!}{k}.$ (15)

The determinant of the matrix $J_{\mathcal{R},k}$ can then be written as:

$\displaystyle\det\quantity[J_{\mathcal{R},k}]$
$\displaystyle=\quantity(-1)^{k+1}N_{\mathcal{R}}\frac{k!}{k}\ell_{\mathcal{R},k}+N_{\mathcal{R}}\sum_{i=1}^{k-1}\quantity(-1)^{i+1}\ell_{\mathcal{R},i}\det\quantity[J_{\mathcal{R},k}^{(i,1)}].$
(16)

Using the coefficients of the characteristic polynomial of the thermal Wilson
line, defined in Eq. (11), we can write the determinant in the left-hand side
of the above equation as:

$\displaystyle(-1)^{-k}k!c_{k}(L_{\mathcal{R}})$
$\displaystyle=\quantity(-1)^{k+1}N_{\mathcal{R}}\frac{k!}{k}\ell_{\mathcal{R},k}+N_{\mathcal{R}}\sum_{i=1}^{k-1}\quantity(-1)^{i+1}\ell_{\mathcal{R},i}\det\quantity[J_{\mathcal{R},k}^{(i,1)}].$
(17)

At this point we still cannot solve the equation for $\ell_{\mathcal{R},k}$
since the coefficients, $c_{k}(L_{\mathcal{R}})$, still strongly depend on
$\ell_{\mathcal{R},k}$. So, we make use of a special property of these
coefficients. If an $N-$dimensional square matrix $M$ is unitary, then the
following holds:

$\displaystyle c_{j}(M)$
$\displaystyle=\quantity(-1)^{N}\det\quantity[M]c_{N-j}(M^{\dagger}),$ (18)

with $j=0,1,2,3,\dots,N-1,N$. Here, $c_{j}(M)$ and $c_{j}(M^{\dagger})$ are
the $j-$coefficients of the characteristic polynomial of the matrices $M$ and
$M^{\dagger}$, respectively. This property is shown in Appendix A. For our
purposes, it will allow us to write $c_{k}(L_{\mathcal{R}})$ in terms of
$c_{N_{\mathcal{R}}-k}(L_{\mathcal{R}}^{\dagger})$. This might not sound as
much progress however, the highest charge Polyakov loop present in the
coefficient $c_{N_{\mathcal{R}}-k}(L_{\mathcal{R}}^{\dagger})$ is
$\overline{\ell}_{\mathcal{R},N_{\mathcal{R}}-k}$. As long as
$N_{\mathcal{R}}-k>0$, all the Polyakov loops present in
$c_{N_{\mathcal{R}}-k}(L_{\mathcal{R}}^{\dagger})$ will have a smaller charge
than $\ell_{\mathcal{R},k}$. Hence, applying the identity written in Eq. (18)
to the thermal Wilson line, $L_{\mathcal{R}}$, we can write Eq. (17) as:

$\displaystyle(-1)^{N_{\mathcal{R}}-k}k!c_{N_{\mathcal{R}}-k}(L_{\mathcal{R}}^{\dagger})$
$\displaystyle=\quantity(-1)^{k+1}N_{\mathcal{R}}\frac{k!}{k}\ell_{\mathcal{R},k}+N_{\mathcal{R}}\sum_{i=1}^{k-1}\quantity(-1)^{i+1}\ell_{\mathcal{R},i}\det\quantity[J_{\mathcal{R},k}^{(i,1)}].$
(19)

Additionally, we have used the fact that the thermal Wilson line is an special
unitary matrix, $\det\quantity[L_{\mathcal{R}}]=1$. Finally, we can solve this
equation for $\ell_{\mathcal{R},k}$ to yield the generating formula to
decompose a non-Cartan-Polyakov loop into combinations of loops with smaller
charges:

$\displaystyle\ell_{\mathcal{R},k}$
$\displaystyle=(-1)^{k+1}\frac{k}{N_{\mathcal{R}}\quantity(N_{\mathcal{R}}-k)!}\det\quantity[\overline{J}_{\mathcal{R},N_{\mathcal{R}}-k}]-\frac{k}{k!}\sum_{i=1}^{k-1}\quantity(-1)^{i-k}\ell_{\mathcal{R},i}\det\quantity[J_{\mathcal{R},k}^{(i,1)}].$
(20)

In the above, we used Eq. (11) to write,
$c_{N_{\mathcal{R}}-k}(L_{\mathcal{R}}^{\dagger})=\nicefrac{{(-1)^{N_{\mathcal{R}}-k}}}{{\quantity(N_{\mathcal{R}}-k)!}}\det\quantity[\overline{J}_{\mathcal{R},N_{\mathcal{R}}-k}]$
and applied the identity, $(-1)^{2N_{\mathcal{R}}-3k-1}=(-1)^{k+1}$. The most
important feature of this equation is that the right-hand side does not depend
on the non-Cartan-Polyakov loop of charge $k$. Indeed, both
$\det\quantity[\overline{J}_{\mathcal{R},N_{\mathcal{R}}-k}]$ and
$\det\quantity[J_{\mathcal{R},k}^{(i,1)}]$, depend only on Polyakov loops with
charge smaller than $k$.

One can also derive a very similar formula using the characteristic polynomial
of $-L_{\mathcal{R}}$, instead of $L_{\mathcal{R}}$ i.e.,
$p\quantity(-L_{\mathcal{R}},\lambda)$. This polynomial can be written as:

$\displaystyle p\quantity(-L_{\mathcal{R}},\lambda)$
$\displaystyle=\det\quantity[-L_{\mathcal{R}}-\lambda I]$
$\displaystyle=\quantity(-1)^{N_{\mathcal{R}}}\quantity(\lambda^{N_{\mathcal{R}}}+\sum_{i=1}^{{N_{\mathcal{R}}}-1}c_{i}(-L_{\mathcal{R}})\lambda^{{N_{\mathcal{R}}}-i}+\quantity(-1)^{N_{\mathcal{R}}}\det\quantity[-L_{\mathcal{R}}]).$
(21)

Similar to the previous characteristic polynomial, the coefficients are given
by:

$\displaystyle c_{k}(-L_{\mathcal{R}})$
$\displaystyle=\frac{(-1)^{k}}{k!}\det\quantity[R_{\mathcal{R},k}].$ (22)

The $k\times k$ matrix $R_{\mathcal{R},k}$ used above, is defined by:

The matrix elements are given by:

$\displaystyle(R_{\mathcal{R},k})_{i,j}=(-1)^{i-j+1}N_{\mathcal{R}}\ell_{\mathcal{R},i-j+1}H\quantity(i-j)+i\delta_{i,j-1}.$
(24)

Following the same steps as before, one can arrive at the following
decomposition formula:

$\displaystyle\ell_{\mathcal{R},k}$
$\displaystyle=(-1)^{N_{\mathcal{R}}+1}\frac{k}{N_{\mathcal{R}}(N_{\mathcal{R}}-k)!}\det\quantity[\overline{R}_{\mathcal{R},N_{\mathcal{R}}-k}]-\frac{k}{k!}\sum_{i=1}^{k-1}\ell_{\mathcal{R},i}\det\quantity[R_{\mathcal{R},k}^{(i,1)}].$
(25)

Once again, the matrix $\overline{R}_{\mathcal{R},k}$, has the same structure
as the matrix defined in Eq. (LABEL:RRepk_matrix_def) with the replacement:
$\ell_{\mathcal{R},k}\to\overline{\ell}_{\mathcal{R},k}$.

We highlight that both formulas, given in Eqs. (20) and (25), are only valid
for $k<N_{\mathcal{R}}$. In order to decompose non-Cartan loops with charge
$k\geq N_{\mathcal{R}}$, one can use the Cayley-Hamilton theorem. Using this
theorem, we can write:

$\displaystyle
L_{\mathcal{R}}^{N_{\mathcal{R}}}+\sum_{i=1}^{{N_{\mathcal{R}}}-1}c_{i}(L_{\mathcal{R}})L_{\mathcal{R}}^{{N_{\mathcal{R}}}-i}+\quantity(-1)^{N_{\mathcal{R}}}=0.$
(26)

Multiplying by $p-$powers of the thermal Wilson line, dividing by the
dimension of the representation and taking the trace, yields the following
equality:

$\displaystyle\frac{1}{N_{\mathcal{R}}}\tr[L_{\mathcal{R}}^{N_{\mathcal{R}}+p}]=-\frac{1}{N_{\mathcal{R}}}\sum_{i=1}^{{N_{\mathcal{R}}}-1}c_{i}(L_{\mathcal{R}})\tr[L_{\mathcal{R}}^{{N_{\mathcal{R}}}+p-i}]+\quantity(-1)^{N_{\mathcal{R}}+1}\frac{1}{N_{\mathcal{R}}}\tr[L_{\mathcal{R}}^{p}].$
(27)

Hence, recognizing the different Polyakov loops, for $k\geq
N_{\mathcal{R}}+p$, we can write:

$\displaystyle\ell_{\mathcal{R},k}=\quantity(-1)^{N_{\mathcal{R}}+1}\ell_{\mathcal{R},k-N_{\mathcal{R}}}-\sum_{i=1}^{{N_{\mathcal{R}}}-1}c_{i}(L_{\mathcal{R}})\ell_{\mathcal{R},k-i}.$
(28)

This non-Cartan-Polyakov loop decomposition equation is valid for $k\geq
N_{\mathcal{R}}$. One can also use the identity given in Eq. (18) in order to
introduce the adjoint Polyakov loops.

In order to show some examples in which Eqs. (20) and (28) are useful, lets
consider Polyakov loops in the fundamental representation of the
$\mathrm{SU}(N_{c})$ group, with increasing number of colors. Consider first
$N_{c}=3$. For this particular case, the Cartan loops are given by the set
$\mathcal{C}_{3}=\quantity{\ell_{\mathrm{F},1},\overline{\ell}_{\mathrm{F},1}}$.
Polyakov loops with higher charge (so-called non-Cartan-Polyakov loops) can be
written as functions of the loops in this set. The charge 2 loop,
$\ell_{\mathrm{F},2}$, can be calculated using Eq. (20). It yields:

$\displaystyle\ell_{\mathrm{F},2}$
$\displaystyle=3\ell_{\mathrm{F},1}^{2}-2\overline{\ell}_{\mathrm{F},1}.$ (29)

For the higher charges, we use Eq. (28):

$\displaystyle\ell_{\mathrm{F},3}$
$\displaystyle=1-\frac{9}{2}\ell_{\mathrm{F},1}^{3}+\frac{9}{2}\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}$
$\displaystyle=1-9\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}+9\ell_{\mathrm{F},1}^{3},$
(30) $\displaystyle\ell_{\mathrm{F},4}$
$\displaystyle=\ell_{\mathrm{F},1}-\frac{9}{2}\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},2}+\frac{3}{2}\ell_{\mathrm{F},2}^{2}+3\ell_{\mathrm{F},1}\ell_{\mathrm{F},3}$
$\displaystyle=4\ell_{\mathrm{F},1}-36\ell_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},1}+6\overline{\ell}_{\mathrm{F},1}^{2}+27\ell_{\mathrm{F},1}^{4},$
(31) $\displaystyle\ell_{\mathrm{F},5}$
$\displaystyle=\ell_{\mathrm{F},2}-\frac{9}{2}\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},3}+\frac{3}{2}\ell_{\mathrm{F},2}\ell_{\mathrm{F},3}+3\ell_{\mathrm{F},1}\ell_{\mathrm{F},4}$
$\displaystyle=15\ell_{\mathrm{F},1}^{2}+81\ell_{\mathrm{F},1}^{5}-135\ell_{\mathrm{F},1}^{3}\overline{\ell}_{\mathrm{F},1}+45\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}^{2}-5\overline{\ell}_{\mathrm{F},1},$
(32) $\displaystyle\;\;\vdots$

Naturally, the non-Cartan-Polyakov loops decomposition formula maintains the
symmetry charge of the loop. For each decomposed non-Cartan-Polyakov loop,
each term in its decomposition, has the same charge as the original non-
Cartan-Polyakov loop. In the case of $\ell_{\mathrm{F},2}$, for example, both
$\ell_{\mathrm{F},1}^{2}$ and $\overline{\ell}_{\mathrm{F},1}$ transform as
fields of charge $+2$. One can readily obtain the respective adjoint Polyakov
loops by replacing, $\ell_{\mathcal{R},k}\to\overline{\ell}_{\mathcal{R},k}$
and $\overline{\ell}_{\mathcal{R},k}\to\ell_{\mathcal{R},k}$.

Increasing the number of colors to $N_{c}=5$ changes the set that constitute
the Cartan loops. In this new case the set is:
$\mathcal{C}_{5}=\quantity{\ell_{\mathrm{F},1},\ell_{\mathrm{F},2},\overline{\ell}_{\mathrm{F},1},\overline{\ell}_{\mathrm{F},2}}$.
Once again, using Eq. (20), we can decompose the non-Cartan loops, with charge
$k<5$, in terms of the Cartan set as follows:

$\displaystyle\ell_{\mathrm{F},3}$
$\displaystyle=\frac{15}{2}\overline{\ell}_{\mathrm{F},1}^{2}-\frac{3}{2}\overline{\ell}_{\mathrm{F},2}-\frac{25}{2}\ell_{\mathrm{F},1}^{3}+\frac{15}{2}\ell_{\mathrm{F},1}\ell_{\mathrm{F},2},$
(33) $\displaystyle\ell_{\mathrm{F},4}$
$\displaystyle=\frac{125}{6}\ell_{\mathrm{F},1}^{4}-4\overline{\ell}_{\mathrm{F},1}-25\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},2}+\frac{20}{3}\ell_{\mathrm{F},1}\ell_{\mathrm{F},3}+\frac{5}{2}\ell_{\mathrm{F},2}^{2}$
$\displaystyle=50\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}^{2}-10\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},2}-4\overline{\ell}_{\mathrm{F},1}-\frac{125}{2}\ell_{1}^{4}+25\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},2}+\frac{5}{2}\ell_{\mathrm{F},2}^{2}.$
(34)

Again, the higher charge Polyakov loops can be decomposed using Eq. (28):

$\displaystyle\ell_{\mathrm{F},5}$
$\displaystyle=1-\frac{625}{24}\ell_{\mathrm{F},1}^{5}+\frac{625}{12}\ell_{\mathrm{F},1}^{3}\ell_{\mathrm{F},2}-\frac{125}{6}\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},3}-\frac{125}{8}\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}^{2}+\frac{25}{4}\ell_{\mathrm{F},1}\ell_{\mathrm{F},4}+\frac{25}{6}\ell_{\mathrm{F},2}\ell_{\mathrm{F},3}$
$\displaystyle=1+\frac{625}{4}\ell_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},1}^{2}-\frac{125}{4}\ell_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},2}-25\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}+\frac{125}{4}\ell_{\mathrm{F},2}\overline{\ell}_{\mathrm{F},1}^{2}-\frac{25}{4}\ell_{\mathrm{F},2}\overline{\ell}_{\mathrm{F},2}-\frac{625}{4}\ell_{\mathrm{F},1}^{5}+\frac{125}{4}\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}^{2},$
(35) $\displaystyle\ell_{\mathrm{F},6}$
$\displaystyle=\ell_{\mathrm{F},1}+5\ell_{\mathrm{F},1}\ell_{\mathrm{F},5}-\frac{625}{24}\ell_{\mathrm{F},2}\ell_{\mathrm{F},1}^{4}+\frac{125}{6}\ell_{\mathrm{F},1}^{3}\ell_{\mathrm{F},3}+\frac{125}{4}\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},2}^{2}$
$\displaystyle\quad\quad\quad\quad\quad\quad-\frac{25}{2}\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},4}-\frac{125}{6}\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}\ell_{\mathrm{F},3}+\frac{5}{3}\ell_{\mathrm{F},3}^{2}+\frac{15}{4}\ell_{\mathrm{F},2}\ell_{\mathrm{F},4}-\frac{25}{8}\ell_{\mathrm{F},2}^{3}$
$\displaystyle=6\ell_{\mathrm{F},1}-75\ell_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},1}+375\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}^{2}\ell_{\mathrm{F},2}-75\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}\overline{\ell}_{\mathrm{F},2}+\frac{375}{4}\overline{\ell}_{\mathrm{F},1}^{4}$
$\displaystyle\quad\quad\quad\quad\quad\quad+\frac{15}{4}\overline{\ell}_{\mathrm{F},2}^{2}-15\overline{\ell}_{\mathrm{F},1}\ell_{\mathrm{F},2}-\frac{75}{2}\overline{\ell}_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},2}-\frac{1875}{4}\ell_{\mathrm{F},2}\ell_{\mathrm{F},1}^{4}+\frac{375}{2}\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},2}^{2}+\frac{25}{4}\ell_{\mathrm{F},2}^{3},$
(36) $\displaystyle\;\;\vdots$

The formalism developed earlier allows one to go on decomposing non-Cartan-
Polyakov loops for even higher values of $N_{c}$. Naturally, as $N_{c}$
increases, the number of terms in each decomposition gets larger and larger.
For $N_{c}=9$, the set is
$\mathcal{C}_{9}=\quantity{\ell_{\mathrm{F},1},\ell_{\mathrm{F},2},\ell_{\mathrm{F},3},\ell_{\mathrm{F},4},\overline{\ell}_{\mathrm{F},1},\overline{\ell}_{\mathrm{F},2},\overline{\ell}_{\mathrm{F},3},\overline{\ell}_{\mathrm{F},4}}$,
for $N_{c}=11$, one must also join $\ell_{\mathrm{F},5}$ and
$\overline{\ell}_{\mathrm{F},5}$ and so on. The number of terms that appear in
the decomposition of non-Cartan-Polyakov loops into Cartan-Polyakov loops
naturally increases with the number of colors. For example, for $N_{c}=11$,
the first non-Cartan loop, $\ell_{\mathrm{F},6}$, has 17 terms while, for
$N_{c}=21$, the first non-Cartan is $\ell_{\mathrm{F},11}$, which contains 97
different terms.

These relations can be checked by hand if one considers a particular gauge
choice. In the so-called Polyakov loop gauge, the $A_{4}$ field in the
fundamental representation is diagonal and static [43] and it can be expressed
by the $N_{c}$ eigenphases $q_{i}$,
$A_{4}=\operatorname{diag}\quantity(iq_{1},iq_{2},\ldots,iq_{N_{c}})$ [44].
Since $A_{4}$ is traceless, the constraint $\sum_{i=1}^{N_{c}}q_{i}=0$
applies. Thus, in this particular gauge, Polyakov loops of different charges
can be written as:

$\displaystyle\ell_{\mathrm{F},k}$
$\displaystyle=\frac{1}{N_{c}}\sum_{j=1}^{N_{c}}\mathrm{e}^{ikq_{j}}=\frac{1}{N_{c}}\quantity(\sum_{j=1}^{N_{c}-1}\mathrm{e}^{ikq_{j}}+\mathrm{e}^{-\sum_{j=1}^{N_{c}-1}ikq_{j}}),$
(37) $\displaystyle\overline{\ell}_{\mathrm{F},k}$
$\displaystyle=\frac{1}{N_{c}}\sum_{j=1}^{N_{c}}\mathrm{e}^{-ikq_{j}}=\frac{1}{N_{c}}\quantity(\sum_{j=1}^{N_{c}-1}\mathrm{e}^{-ikq_{j}}+\mathrm{e}^{\sum_{j=1}^{N_{c}-1}ikq_{j}}).$
(38)

Using these gauge fixed Polyakov loops in the fundamental representation, one
can check the results coming from the non-Cartan-Polyakov decomposition
formulas.

##  IV Applications

In this section we show two examples of how to construct models of quarks and
gluons with any number of colors. First, we consider the PNJL model in the
mean field approximation. Secondly, we consider a glue effective potential.

###  IV.1 The mean-field Polyakov$-$Nambu$-$Jona-Lasinio model with an
arbitrary $N_{c}$

The NJL model is widely used as an effective model of QCD due to its
simplicity, symmetry properties and ability to display dynamical chiral
symmetry breaking in the vacuum and restoration at high
temperatures/densities. Some applications of this model include the study of
in-medium meson behavior, transport coefficients of quark matter, the equation
of state of hybrid neutron stars and the phase diagram of strongly interacting
matter with great emphasis on the location of the conjectured critical
endpoint. This model however, lacks the ability to describe the confinement-
deconfinement transition. In the Polyakov version of the model, a static
gluonic background field is coupled to the quarks sector, allowing for the
description of both the chiral and deconfinement transitions [33, 45, 46, 47].

The Polyakov loop can be included in the model by minimally coupling a
background gluon field in the time direction, $A_{0}=-iA_{4}$ [46], to the
Lagrangian density via the covariant derivative. It can be defined as:

$\displaystyle\mathcal{L}\quantity(\psi,\overline{\psi},A_{4},T)$
$\displaystyle=\overline{\psi}\quantity(i\not{D}-\hat{m})\psi+\mathcal{L}_{\mathrm{int}}\quantity(\psi,\overline{\psi})-\mathcal{U}_{\mathrm{eff}}\quantity(A_{4},T).$
(39)

Here, $\psi$ is the quark field,
$D^{\mu}=\partial^{\mu}-i\delta_{0}^{\mu}A^{0}$ is the covariant derivative
and $\hat{m}=\operatorname{diag}\quantity{m_{1},m_{2},\ldots,m_{N_{f}}}$ is
the quark current mass matrix (with $N_{f}$ the number of quark flavors). The
different quark-quark interactions are contained in the term
$\mathcal{L}_{\mathrm{int}}\quantity(\psi,\overline{\psi})$ and, for our
purposes, its exact definite form is not important. This interaction term not
only includes the dynamical chiral symmetry breaking 4-quark scalar-
pseudoscalar interaction,
$\mathcal{L}\supset(\overline{\psi}\lambda_{a}\psi)^{2}+(\overline{\psi}i\gamma_{5}\lambda_{a}\psi)^{2}$
($\lambda_{a}$ are the generators of the $\mathrm{U}(N_{f})$ algebra), but it
also includes other multi-quark interactions, like the ’t Hooft determinant,
eight quark-quark interactions, explicit chiral symmetry breaking
interactions, vector interactions, etc [48, 49, 50, 51]. The effective
potential, $\mathcal{U}_{\mathrm{eff}}\quantity(A_{4},T)$, incorporates the
spontaneous breaking of the $\mathrm{Z}(N_{c})$ symmetry at some finite
temperature. This potential will be discussed in much detail in later sections
where we will define a version of it for any number of colors.

The Euclidean generating functional for the PNJL model, can be written as:

$\displaystyle\mathcal{Z}\quantity[\eta,\overline{\eta},j_{4}]\propto\int\mathcal{D}A_{4}\mathcal{D}\overline{\psi}\mathcal{D}\psi\exp\quantity[-\mathcal{S}_{E}\quantity[\overline{\psi},\psi,A_{4}]+\overline{\psi}\eta+\overline{\eta}\psi+A_{4}j_{4}].$
(40)

Here, $\mathcal{S}_{E}$ is the Euclidean action of the PNJL model. The path
integral over the temporal gluon field $A_{4}$ is to be made over all of its
$N_{c}-1$ components and $j_{4}$ is a source for the temporal gluon field.

To deal with the path integral over the quark fields, we apply the mean field
approximation. In this scheme, the multi-quark interactions present in
$\mathcal{L}_{\mathrm{int}}$, can be written as the product of quark bilinear
operators with the introduction of mean fields, $\expectationvalue{\phi_{i}}$.
After the linearization of the product between these operators, a quadratic
Lagrangian in the quark fields is obtained and the quark fields can be
integrated out from the generating function [49]. The generating functional
becomes:

$\displaystyle\mathcal{Z}\quantity[\expectationvalue{\phi_{i}},j_{4}]\propto\int\mathcal{D}A_{4}\exp\quantity[-\mathcal{S}_{E}\quantity[\expectationvalue{\phi_{i}},A_{4}]+A_{4}j_{4}].$
(41)

The only thing left to do is to solve the path integral over the temporal
gluon field. Different approaches to deal with the path integral over the
$A_{4}$ field, can be found in the literature, ranging from using the
Functional Renormalization Group, a Gaussian model or using the Weiss mean
field approximation, see [41, 52, 53]. One can parametrize the
$\mathrm{SU}(N_{c})$ group volume in terms of the diagonal elements of the
$\mathrm{SU}(N_{c})$ algebra, i.e., the Cartan elements of the algebra [54,
55, 56, 57]. Hence, we can separate the integration over $A_{4}$ in Eq. (41)
in two contributions, one made over the diagonal generators of the
$\mathrm{SU}(N_{c})$ algebra and another over the non-diagonal generators. One
can integrate over the non-diagonal components to yield the so-called Haar
measure of the $\mathrm{SU}(N_{c})$ group, $H$ [58, 55, 59]:

$\displaystyle
H(A_{4}^{\operatorname{diag}})=\frac{1}{V_{\mathrm{SU}(N_{c})}}\int\prod_{i\notin\operatorname{diag}}\differential{A_{4}^{(i)}}.$
(42)

Here, $A_{4}^{\mathrm{diag}}$ is the set of diagonal components of $A_{4}$. In
the already introduced Polyakov gauge, $A_{4}$ is diagonal and traceless,
meaning it can be parametrized in terms of the Cartan elements,
$A_{4}=\sum_{i\in\operatorname{diag}}A_{4}^{(i)}t_{i}$ with $t_{i}$ the
diagonal generators of the $\mathrm{SU}(N_{c})$ algebra i.e., the Cartan
subalgebra. For $N_{c}=3$, for instance, we can write
$A_{4}=A_{4}^{(3)}\lambda_{3}+A_{4}^{(8)}\lambda_{8}$, with $\lambda_{3}$ and
$\lambda_{8}$ the diagonal Gell-Mann matrices. Thus, in Eq. (41) one can
integrate out the non-diagonal elements of $A_{4}$ leaving only the path
integral over the diagonal components, which contains all the physics of the
model [58, 55]:

$\displaystyle\mathcal{Z}\quantity[\expectationvalue{\phi_{i}},j_{4}]\propto\int\prod_{j\in\operatorname{diag}}\mathcal{D}{A_{4}^{(j)}}\exp\quantity[-\mathcal{S}_{E}\quantity[\expectationvalue{\phi_{i}},A_{4}]+A_{4}j_{4}+\kappa\ln\quantity[H(A_{4}^{\operatorname{diag}})]].$
(43)

Here, $\kappa$ is some proportionality constant and
$H[A_{4}^{\operatorname{diag}}]$ is the Haar measure, which results from
integrating the non-diagonal components of $A_{4}$. The Haar measure can be
absorbed on the definition of the Polyakov effective potential,
$\mathcal{U}_{\mathrm{eff}}$. As a matter of fact, the Haar measure has been
understood to be a crucial part of the potential by favoring the confined
phase [34].

At this point one needs to solve the path integral over the diagonal
components of the field $A_{4}$. The standard approach is to consider that the
functional integration is dominated by the stationary point: quantum
fluctuations are neglected and only the classical configuration contributes to
the path integral [60, 50]. Thus, the functional integration can be dropped
and the diagonal components of the $A_{4}$ field can be treated as mean
fields. The values of these fields are calculated by requiring the effective
action to be stationary, $\derivative*{\mathcal{S}_{E}}{A_{4}^{(i)}}=0$, with
$A_{4}^{(i)}$ the Cartan elements.

Usually, for $N_{c}=3$, the requirements of stationarity of the effective
action with respect to the diagonal components of the $A_{4}$ field are
replaced by stationarity with respect to the values of the charge one Polyakov
loop and its adjoint in the fundamental representation, $\Phi$ and
$\overline{\Phi}$ [60]. Since these quantities are functions of
$A_{4}^{\operatorname{diag}}$, $\Phi=\Phi[A_{4}^{\operatorname{diag}}]$ and
$\overline{\Phi}=\overline{\Phi}[A_{4}^{\operatorname{diag}}]$, changing
variables is allowed as long as the Jacobian of the transformation is non-
singular [61, 62, 63]. For other values of $N_{c}$, $\ell_{\mathrm{F},1}=\Phi$
and $\overline{\ell}_{\mathrm{F},1}=\overline{\Phi}$ are not enough degrees of
freedom to express the Haar measure. Thus, for other values of $N_{c}$, we
will consider that these change of variables are allowed and replace the
variables of the model: instead of using the set of Cartan elements of the
$A_{4}$ field, we will use the Cartan-Polyakov loop set defined in Eq. (9).

Within the mean field approximation, the Lagrangian density of the PNJL model
can be written as the one of fermionic quasiparticles with an effective mass,
$M$, chemical potential, $\mu$, and effective mean field potential, $U$. Using
the Matsubara formalism, one can calculate the thermodynamic potential of the
model at finite temperature and chemical potential,
$\Omega\quantity(N_{c},T,\mu)$ [33, 46]. It is given by:

$\displaystyle\Omega(N_{c},T,\mu)$
$\displaystyle=U(\expectationvalue{\phi_{i}})+\mathcal{U}_{\mathrm{eff}}(A_{4}^{\operatorname{diag}},T)$
$\displaystyle-2T\tr\int_{\mathrm{reg}}\frac{\differential[3]{p}}{\quantity(2\pi)^{3}}\quantity{\beta
N_{c}E+\ln\det\quantity[1+L_{\mathrm{F}}\mathrm{e}^{-\beta\quantity(E-\mu)}]+\ln\det\quantity[1+L_{\mathrm{F}}^{\dagger}\mathrm{e}^{-\beta\quantity(E+\mu)}]}.$
(44)

Here, $E=\sqrt{p^{2}+M^{2}}$ and $\beta$ is the inverse temperature,
$\beta={T}^{-1}$. The potential $U(\expectationvalue{\phi_{i}})$ contains the
contribution coming from the mean fields, $\expectationvalue{\phi_{i}}$. The
trace must be calculated over flavor indexes while the determinant is
calculated over color space. The mean fields are fixed by requiring for
thermodynamic consistency,
$\derivative*{\Omega}{\expectationvalue{\phi_{i}}}=0$ [64]. The subscript
“reg” in the integration, serves only as an acknowledgment for the fact that
the PNJL model has divergent integrals and some regularization scheme must be
employed to render these integrals finite. Some techniques include the
3-momentum or the Pauli-Villars regularization schemes.

Until now, we have not yet specified the number of colors. In order to obtain
the thermodynamic potential for a particular number of colors, one must write
the effective Polyakov loop potential and calculate the determinant inside the
logarithms, for the chosen $N_{c}$.

####  IV.1.1 The fermionic determinant with an arbitrary $N_{c}$

Lets start by evaluating the fermionic contribution for a given number of
colors. The fermionic determinants are given by:

$\displaystyle\det\quantity[1+L_{\mathrm{F}}\mathrm{e}^{-\beta\quantity(E-\mu)}]=h_{-}^{N_{c}}\det\quantity[h_{-}^{-1}I+L_{\mathrm{F}}],$
(45)
$\displaystyle\det\quantity[1+L_{\mathrm{F}}^{\dagger}\mathrm{e}^{-\beta\quantity(E+\mu)}]=h_{+}^{N_{c}}\det\quantity[h_{+}^{-1}I+L_{\mathrm{F}}^{\dagger}].$
(46)

Here, we have defined the non-zero quantities:

$\displaystyle h_{\pm}=\mathrm{e}^{-\beta\quantity(E\pm{\mu})}.$ (47)

The goal is to calculate the determinants for an arbitrary $N_{c}$ using only
the corresponding set of Cartan-Polyakov loops. To this end, we will write the
determinants in Eqs. (45) and (46) in their respective characteristic
polynomials. These polynomials will have $(N_{c}-1)$ coefficients,
$c_{k}(-M)$, that have to be calculated. Calculating the $k-$th coefficient,
implies the calculation of another determinant, as laid out in Eq. (11).
Considering the characteristic polynomial
$p\quantity(-L_{\mathcal{R}},h_{-}^{-1})$ (see Eq. (21)), we can write Eq.
(45) as:

$\displaystyle h_{-}^{N_{c}}\det\quantity[L_{\mathrm{F}}+h_{-}^{-1}I]$
$\displaystyle=1+\sum_{k=1}^{{N_{c}}-1}c_{k}(-L_{\mathrm{F}})h_{-}^{k}+h_{-}^{N_{c}}.$
(48)

We used the fact the dimension of the fundamental representation is,
$N_{\mathrm{F}}=N_{c}$ and that
$\det\quantity[-L_{\mathrm{F}}]=\quantity(-1)^{N_{c}}$.

In the sum, the difficulty of calculating the determinant increases with the
number of colors. For a fixed number of colors, calculating the $k-$th
coefficient in the sum, amounts to solving the determinant of the
$R_{\mathcal{R},k}$ matrix, see Eq. (22). Not only that, the $k-$th
coefficient in the sum will necessarily contain Polyakov loops, of charge
$\quantity[1,\ldots,k]$. This can be seen by observing the structure of the
matrix $R_{\mathcal{R},k}$, defined in Eq. (LABEL:RRepk_matrix_def). As
introduced in this work, by definition, the Cartan-Polyakov loops have charge
in the interval $\quantity[1,\ldots,(N_{c}-1)/2]$. Hence, the coefficients in
the sum, with $k=(N_{c}-1)/2+1,\ldots,N_{c}-1$, will unavoidably contain non-
Cartan-Polyakov loops. Since the goal is to write the model exclusively as a
function of Cartan-Polyakov loops, this feature poses the difficult task of
expanding the determinant of Eq. (48) and then applying the charge
decomposition formula, given in Eq. (20), to systematically decompose all non-
Cartan-Polyakov loops into Cartan ones. In order to solve this issue, we will
use a trick to simplify the practical calculations while, at the same time,
automatically re-write the sum as a function of Cartan-Polyakov loops only.
For simplicity, in the following discussion, we will focus on the case in
which the number of colors is odd. The even number of colors case is a
straightforward extension.

Consider the identity given in Eq. (18), which states that we can relate the
coefficient $c_{k}(-L_{\mathrm{F}})$, with the coefficient
$c_{N_{c}-k}(-L_{\mathrm{F}}^{\dagger})$. If we employ this identity, we can
break the original sum in Eq. (48) in the middle, in two distinct sums and
write,

$\displaystyle\sum_{k=1}^{{N_{c}}-1}c_{k}(-L_{\mathrm{F}})h_{-}^{k}=\sum_{k=1}^{(N_{c}-1)/2}c_{k}(-L_{\mathrm{F}})h_{-}^{k}+\sum_{k=(N_{c}+1)/2}^{{N_{c}}-1}c_{N_{c}-k}(-L_{\mathrm{F}}^{\dagger})h_{-}^{k},$
(49)

where we have used
$c_{j}(-L_{\mathrm{F}})=c_{N-j}(-L_{\mathrm{F}}^{\dagger})$. One can check
that both the first and second sums, on the right-hand side of this equation,
contains coefficients that will only produce Polyakov loops with charges in
the interval $\quantity[1,\ldots,(N_{c}-1)/2]$. This range of charges is the
same as the one which defines the Cartan-Polyakov loop set. Hence, the
determinant can be written only as a function of Cartan-Polyakov loops.

From an analytically point of view, another nice feature of this trick is that
by employing it, one only needs to calculate half of the coefficients. At
first glance, one has to calculate $N_{c}-1$ coefficients, starting with the
coefficient with $k=1$ up to $k=N_{c}-1$. However, one can calculate the
coefficients, $c_{k}(-L_{\mathrm{F}})$, up to $k=(N_{c}-1)/2$ and the
remaining coefficients, related to $c_{k}(-L_{\mathrm{F}}^{\dagger})$, can be
obtained from the substitution,
$\ell_{\mathrm{F},k}\to\overline{\ell}_{\mathrm{F},k}$.

Thus, for an odd number of colors, we can re-write the sum using Eq. (49), and
write the determinant in Eq. (45) as:

$\displaystyle\det\quantity[1+L_{\mathrm{F}}\mathrm{e}^{-\beta\quantity(E-\mu)}]$
$\displaystyle=1+\sum_{k=1}^{(N_{c}-1)/2}c_{k}(-L_{\mathrm{F}})h_{-}^{k}+\sum_{k=(N_{c}+1)/2}^{{N_{c}}-1}c_{N_{c}-k}(-L_{\mathrm{F}}^{\dagger})h_{-}^{k}+h_{-}^{N_{c}}.$
(50)

For the case in which $N_{c}$ is an even number, the idea is the same with the
small caveat that one of the sums will have an additional term.

The calculation of the determinant, defined in Eq. (46), follows the same
steps. One gets the same result as the one shown in Eq. (50) with the
following replacements: $-\to+$, $L_{\mathrm{F}}\to L_{\mathrm{F}}^{\dagger}$
and $L_{\mathrm{F}}^{\dagger}\to L_{\mathrm{F}}$.

The explicit calculation of the coefficients follows from Eq. (22). Hence, one
can readily calculate the determinants defined in Eqs. (45) and (46), for an
arbitrary $N_{c}$. The problem of writing the fermionic part of the PNJL model
for an arbitrary number of colors is reduced to calculating the
$c_{k}(-L_{\mathrm{F}})$ coefficients up to $(N_{c}-1)/2$. The first
coefficients are:

$\displaystyle c_{1}(-L_{\mathrm{F}})$
$\displaystyle=N_{c}\ell_{\mathrm{F},1},$ (51) $\displaystyle
c_{2}(-L_{\mathrm{F}})$
$\displaystyle=\frac{N_{c}}{2}\quantity(N_{c}\ell_{\mathrm{F},1}^{2}-\ell_{\mathrm{F},2}),$
(52) $\displaystyle c_{3}(-L_{\mathrm{F}})$
$\displaystyle=\frac{N_{c}}{6}\quantity(N_{c}^{2}\ell_{\mathrm{F},1}^{3}-3N_{c}\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}+2\ell_{\mathrm{F},3}),$
(53) $\displaystyle c_{4}(-L_{\mathrm{F}})$
$\displaystyle=\frac{N_{c}}{24}\quantity(N_{c}^{3}\ell_{\mathrm{F},1}^{4}-6N_{c}^{2}\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},2}+3N_{c}\ell_{\mathrm{F},2}^{2}+8N_{c}\ell_{\mathrm{F},1}\ell_{\mathrm{F},3}-6\ell_{\mathrm{F},4}),$
(54) $\displaystyle\;\;\vdots$

As an example, for $N_{c}=9$, the fermionic determinant containing
$L_{\mathrm{F}}$, is given by:

$\displaystyle\det\quantity[1+L_{\mathrm{F}}\mathrm{e}^{-\beta\quantity(E-\mu)}]$
$\displaystyle=1$
$\displaystyle+9\ell_{\mathrm{F},1}\mathrm{e}^{-\beta\quantity(E-\mu)}$
$\displaystyle+\frac{9}{2}\quantity(9\ell_{\mathrm{F},1}^{2}-\ell_{\mathrm{F},2})\mathrm{e}^{-2\beta\quantity(E-\mu)}$
$\displaystyle+\frac{3}{2}\quantity(81\ell_{\mathrm{F},1}^{3}-27\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}+2\ell_{\mathrm{F},3})\mathrm{e}^{-3\beta\quantity(E-\mu)}$
$\displaystyle+\frac{9}{8}\quantity(243\ell_{\mathrm{F},1}^{4}-162\ell_{\mathrm{F},1}^{2}\ell_{\mathrm{F},2}+9\ell_{\mathrm{F},2}^{2}+24\ell_{\mathrm{F},1}\ell_{\mathrm{F},3}-2\ell_{\mathrm{F},4})\mathrm{e}^{-4\beta\quantity(E-\mu)}$
$\displaystyle+\frac{9}{8}\quantity(243\overline{\ell}_{\mathrm{F},1}^{4}-162\overline{\ell}_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},2}+9\overline{\ell}_{\mathrm{F},2}^{2}+24\overline{\ell}_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},3}-2\overline{\ell}_{\mathrm{F},4})\mathrm{e}^{-5\beta\quantity(E-\mu)}$
$\displaystyle+\frac{3}{2}\quantity(81\overline{\ell}_{\mathrm{F},1}^{3}-27\overline{\ell}_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},2}+2\overline{\ell}_{\mathrm{F},3})\mathrm{e}^{-6\beta\quantity(E-\mu)}$
$\displaystyle+\frac{9}{2}\quantity(9\overline{\ell}_{\mathrm{F},1}^{2}-\overline{\ell}_{\mathrm{F},2})\mathrm{e}^{-7\beta\quantity(E-\mu)}$
$\displaystyle+9\overline{\ell}_{\mathrm{F},1}\mathrm{e}^{-8\beta\quantity(E-\mu)}$
$\displaystyle+\mathrm{e}^{-9\beta\quantity(E-\mu)}.$ (55)

The determinant containing $L_{\mathrm{F}}^{\dagger}$, can be obtained from
the above expression with the appropriate substitutions:
$\ell_{\mathrm{F},k}\to\overline{\ell}_{\mathrm{F},k}$,
$\overline{\ell}_{\mathrm{F},k}\to\ell_{\mathrm{F},k}$ and $\mu\to-\mu$.

####  IV.1.2  The Polyakov loop effective potential with an arbitrary $N_{c}$

We now turn our attention to the effective Polyakov loop potential,
$\mathcal{U}_{\mathrm{eff}}$, which can be written using the Ginzburg$-$Landau
theory of phase transitions. Within this approach, the potential must respect
the symmetry properties of the original theory.

For $N_{c}=3$, one can build this potential using the charge one Polyakov loop
in the fundamental representation (and its conjugate),
$\ell_{\mathrm{F},1}=\Phi$ and
$\overline{\ell}_{\mathrm{F},1}=\overline{\Phi}$. Different phenomenological
potentials were introduced in the literature including the logarithmic,
polynomial and polynomial-logarithmic. They all share a kinetic term,
proportional to $\Phi\overline{\Phi}$ and higher order terms that are
essential to build a bounded potential, with the correct symmetries (the
kinetic term introduces an additional $\mathrm{U}(1)$ symmetry which must be
removed by other terms [54, 56, 65]), which reproduces a first-order phase
transition at some finite temperature.

In the logarithm potential, the non-quadratic term is given by the logarithm
of the Haar measure. As already seen, this contribution naturally appears when
one integrates out the non-diagonal components of the $A_{4}$ field [33, 58].
This potential, for $N_{c}=3$, is widely defined as [33, 60, 66]:

$\displaystyle\frac{\mathcal{U}_{\mathrm{eff}}\quantity(\Phi,\overline{\Phi},T)}{T^{4}}=-\frac{1}{2}a\quantity(T)\overline{\Phi}\Phi+b\quantity(T)\ln\quantity[1-6\overline{\Phi}\Phi+4\quantity(\overline{\Phi}^{3}+\Phi^{3})-3\quantity(\overline{\Phi}\Phi)^{2}].$
(56)

The quantity inside the logarithm is the Haar measure for $N_{c}=3$. The
temperature dependent parameters [60, 66] are:

$\displaystyle a\quantity(T)$
$\displaystyle=a_{0}+a_{1}\quantity(\frac{T_{0}}{T})+a_{2}\quantity(\frac{T_{0}}{T})^{2},$
(57) $\displaystyle b\quantity(T)$
$\displaystyle=b_{3}\quantity(\frac{T_{0}}{T})^{3}.$ (58)

The parameters $T_{0}$, $a_{0}$, $a_{1}$, $a_{2}$ and $a_{3}$ are fixed by
reproducing lattice QCD results [67, 68] with $T_{0}$ fixing the temperature
of the phase transition [46]. The parameter $a_{0}$ is fixed by requiring the
Steffan-Boltzmann limit when $T\to\infty$ while the parameter $b_{3}$ is fixed
by requiring the second-order phase transition to occur at $T=T_{0}$ [46, 60].
The gluon contribution to the Steffan-Boltzmann pressure is
$\nicefrac{{P_{\mathrm{SB}}}}{{T^{4}}}=\nicefrac{{2(N_{c}^{2}-1)\pi^{2}}}{{90}}$.
Hence, the $a_{0}$ parameter can be readily defined at diverging temperatures,
$T\to\infty$, $a_{0}=\nicefrac{{4(N_{c}^{2}-1)\pi^{2}}}{{90}}$. Thus, one can
conclude that this parameter is $N_{c}$ dependent. Since $b_{3}$ can be
calculated from the remaining parameters, including $a_{0}$, it is also
$N_{c}$ dependent. When building an effective Polyakov loop potential for any
number of colors, this has to be taken into account.

For $N_{c}=3$, following lattice QCD calculations at zero chemical potential,
the $T_{0}$ parameter is usually fixed to $270$ MeV [69, 70]. Different
approaches have been used in the literature to fix the value of $T_{0}$,
specially when studying systems at finite chemical potential. As a matter or
fact, in Ref. [71], the parameter $T_{0}$ was given an explicit dependence on
the number of flavors of quarks considered, due to renormalization group
arguments.

The polynomial potential is very similar with the major difference that higher
powers of $\ell_{\mathrm{F},1}=\Phi$ and
$\overline{\ell}_{\mathrm{F},1}=\overline{\Phi}$ are incorporated in the
potential [72, 46]. The polynomial-logarithmic potential, as stated in the
name, is a mixture between the logarithmic and polynomial potentials [73].
Another important property of the logarithmic potential is its diverging
nature in the limit $\ell_{\mathrm{F},1}=\Phi\to 1$ and
$\overline{\ell}_{\mathrm{F},1}=\overline{\Phi}\to 1$. This ensures that the
Polyakov loops do not exceed unity [60].

The question that now arises is the following: how to define an effective
Polyakov loop potential for any given value of the number of colors? The
studies performed in Refs. [26, 28], built effective Polyakov loop potentials
for any number of colors based on the polynomial potential in order to study
the phase diagram at large $N_{c}$. However, such potentials were written only
as a function of $\ell_{\mathrm{F},1}=\Phi$ and
$\overline{\ell}_{\mathrm{F},1}=\overline{\Phi}$ and Polyakov loops with
larger charges were neglected. At large $N_{c}$, the fundamental loop with
charge one may indeed be the most important degree of freedom however, in this
work, we want to write a general expression for any $N_{c}$, without any a
priori approximations (apart from the mean field one). After building such
model, then one may perform a large $N_{c}$ approximation.

In this work, we propose to simply extend the logarithmic potential to any
value of $N_{c}$. Using the polynomial potential, would mean having to
consider different interaction terms between different Cartan-Polyakov loops
in such a way to build $\mathrm{Z}(N_{c})$ invariant and $\mathrm{U}(1)$
breaking interactions. The Haar measure automatically gives interaction terms
which respect the symmetries of the model, yielding an appropriate effective
potential at any number of $N_{c}$. To this end we consider that the potential
can be written solely in terms of the previously introduced Cartan-Polyakov
loops. Hence, for the kinetic part of the potential we simply consider a sum
over quadratic terms of the Cartan-Polyakov loop set:

$\displaystyle\sum_{i=1}^{\frac{N_{c}-1}{2}}\ell_{\mathrm{F},i}\overline{\ell}_{\mathrm{F},i}.$
(59)

This term is invariant under $\mathrm{Z}(N_{c})$ and under an unwanted
$\mathrm{U}(1)$ symmetry which must be removed by the remaining contributions.
The non-quadratic part of the potential will be given by the Haar measure, for
a given $N_{c}$, written exclusively in terms of Cartan-Polyakov loops. In a
later discussion, we will cover the details of how one can accomplish this
task. Hence, the effective Polyakov loop potential, with an arbitrary $N_{c}$,
can be expressed as:

$\displaystyle\frac{\mathcal{U}_{\mathrm{eff}}\quantity(\mathcal{C}_{N_{c}},T)}{T^{4}}=-\frac{1}{2}A\quantity(N_{c},T)\sum_{i=1}^{\frac{N_{c}-1}{2}}\ell_{\mathrm{F},i}\overline{\ell}_{\mathrm{F},i}+B\quantity(N_{c},T)\ln\quantity[H(\mathcal{C}_{N_{c}})].$
(60)

Here, we define the temperature and color dependent parameters
$A\quantity(N_{c},T)$ and $B\quantity(N_{c},T)$, in the same way as it was
done for the temperature parameters $a(T)$ and $b(T)$ in Eqs. (57) and (58):

$\displaystyle A\quantity(N_{c},T)$
$\displaystyle=A_{0}(N_{c})+A_{1}(N_{c})\quantity[\frac{T_{0}}{T}]+A_{2}(N_{c})\quantity[\frac{T_{0}}{T}]^{2},$
(61) $\displaystyle B\quantity(N_{c},T)$
$\displaystyle=B_{3}(N_{c})\quantity[\frac{T_{0}}{T}]^{3}.$ (62)

As before, the $A_{0}(N_{c})$ parameter have to be fixed in such a way to
recover the Steffan-Boltzmann limit for a gas of gluons with
$\mathrm{SU}(N_{c})$ symmetry. The remaining parameters must ensure that a
second-order phase transition occurs at $T_{0}$ at least for $N_{c}=3$. At
higher values of $N_{c}$ other behaviors might be expected [9]. We did not
introduce a color dependence in the $T_{0}$ parameter. Some lattice results,
seem to indicate that $T_{0}$ very weakly depends on the number of colors [74,
75, 76, 9]. Indeed, in such studies, the deconfinement temperature $T_{0}$, as
a function of the number of colors, was found to be,

$\displaystyle T_{0}(N_{c})$
$\displaystyle=\sqrt{\sigma}\quantity(\alpha_{1}+\frac{\alpha_{2}}{N_{c}^{2}}),$
(63)

where, $\sigma$ is the $N_{c}-$independent string tension. In [74, 75], the
values obtained were $\alpha_{1}=0.582(15)$ and $\alpha_{2}=0.43(13)$, in Ref.
[76], $\alpha_{1}=0.9026(23)$ and $\alpha_{2}=0.880(43)$ while, in Ref. [77],
$\alpha_{1}=0.5949(17)$ and $\alpha_{2}=0.458(18)$.

The exact functional dependencies, on the number of colors, of the different
parameters present in the effective Polyakov loop potential, is beyond the
scope of the present work. In order to find such dependencies, one can perform
phenomenological studies involving results coming from lattice QCD with
$N_{c}\geq 2$. The current status of the studies seems to agree that the
transition of pure Yang-Mills theory is of second order for both $N_{c}=2$ and
$N_{c}=3$, while it is of first-order for $N_{c}\geq 4$ [9]. The most
straightforward way to give some $N_{c}$ dependence to the parameters would be
to simply consider the Steffan-Boltzmann requirement enforced on $A_{0}$ and
fix the remaining parameters at $N_{c}=3$.

We now turn our attention to evaluating the argument inside the logarithm of
Eq. (60), $H(\mathcal{C}_{N_{c}})$, for a given $N_{c}$. Since we wrote the
fermionic determinant and the kinetic part of the potential, using Cartan-
Polyakov loops as degrees of freedom, we also want to write the non-quadratic
part of the potential in terms of such variables. In the following, we show
how one can write the argument of the logarithm in Eq. (60), the Haar measure,
in terms of Cartan-Polyakov loops. The Haar measure is given by [78, 34, 79,
44, 59]:

$\displaystyle H(q_{1},q_{2},\ldots,q_{N_{c}})=\prod_{1\leq
i<j}^{N_{c}}\absolutevalue{\mathrm{e}^{iq_{i}}-\mathrm{e}^{iq_{j}}}^{2}.$ (64)

Here, $q_{i}$ are $N_{c}$ eigenphases [44]. This measure can be written using
the Vandermonde determinant. The Vandermonde matrix, $V$, of a set of
variables $z_{1},\ldots,z_{m}$, is defined as [79]:

The matrix elements of this quantity can readily be written:

$\displaystyle\quantity(V)_{ij}$ $\displaystyle=z_{i}^{j-1}.$ (65)

The determinant of this matrix, known as the Vandermonde determinant, can be
calculated to yield:

$\displaystyle\det\quantity[V(z_{1},\ldots,z_{N_{c}})]=\prod_{1\leq i<j\leq
N_{c}}\quantity(z_{j}-z_{i}).$ (66)

The relation between the Vandermonde determinant and the Haar measure is333 We
used the following property:
$\displaystyle\absolutevalue{\det[V\big{]}}^{2}=\det[V^{*}\big{]}\det[V\big{]}=\det[V^{\dagger}\big{]}\det[V\big{]}=\det[V^{\dagger}V\big{]}.$
:

$\displaystyle H(q_{1},q_{2},\ldots,q_{N_{c}})$
$\displaystyle=\absolutevalue{\det[V(\mathrm{e}^{iq_{1}},\ldots,\mathrm{e}^{iq_{N_{c}}})\big{]}}^{2}$
$\displaystyle=\det[V^{\dagger}(\mathrm{e}^{iq_{1}},\ldots,\mathrm{e}^{iq_{N_{c}}})V(\mathrm{e}^{iq_{1}},\ldots,\mathrm{e}^{iq_{N_{c}}})\big{]}$
$\displaystyle=\det[W(q_{1},\ldots,q_{N_{c}})\big{]}.$ (67)

Where,
$W(q_{1},\ldots,q_{N_{c}})=V^{\dagger}(\mathrm{e}^{iq_{1}},\ldots,\mathrm{e}^{iq_{N_{c}}})V(\mathrm{e}^{iq_{1}},\ldots,\mathrm{e}^{iq_{N_{c}}})$.
The matrix elements of $W(q_{1},\ldots,q_{N_{c}})$ are given by:

$\displaystyle\quantity(W)_{ij}$
$\displaystyle=\sum_{l=1}^{N_{c}}\quantity(V^{\dagger})_{il}\quantity(V)_{lj}$
$\displaystyle=\sum_{l=1}^{N_{c}}\quantity(V^{*})_{li}\quantity(V)_{lj}$
$\displaystyle=\sum_{l=1}^{N_{c}}{z^{*}_{l}}^{i-1}z_{l}^{j-1}$
$\displaystyle=\sum_{l=1}^{N_{c}}z_{l}^{j-i}.$ (68)

Here, $z_{l}=\mathrm{e}^{iq_{l}}$. Recognizing that the diagonal entries are
equal to $N_{c}$ and that the non-diagonal entries of this matrix are the
fundamental Polyakov loops with different charges (see Eqs. (37) and (38)),
one can write the matrix $W$ as [79]:

Thus, the Haar measure is exclusively written in terms of Polyakov loops of
different charges (and its adjoints), in the fundamental representation. It is
given by:

$\displaystyle
H(\ell_{\mathrm{F},1},\ldots,\ell_{\mathrm{F},N_{c}-1},\overline{\ell}_{\mathrm{F},1},\ldots,\overline{\ell}_{\mathrm{F},N_{c}-1})$
$\displaystyle=\det\quantity[W].$ (70)

Using the recursion formula in Eq. (20), it is possible to write the Haar
measure above using only the set of Cartan-Polyakov loops,
$\mathcal{C}_{N_{c}}$. Hence,
$H(\ell_{\mathrm{F},1},\ldots,\ell_{\mathrm{F},N_{c}-1},\overline{\ell}_{\mathrm{F},1},\ldots,\overline{\ell}_{\mathrm{F},N_{c}-1})=H(\mathcal{C}_{N_{c}})$.
For $N_{c}=3$, for example, we can write the Haar measure as:

Using Eq. (29), we can write
$\ell_{\mathrm{F},2}=\ell_{\mathrm{F},2}(\ell_{\mathrm{F},1},\overline{\ell}_{\mathrm{F},1})$,
$\overline{\ell}_{\mathrm{F},2}=\overline{\ell}_{\mathrm{F},2}(\ell_{\mathrm{F},1},\overline{\ell}_{\mathrm{F},1})$
and, after the evaluation of the determinant, obtain the same expression as
the one inside the logarithm in Eq. (56). For our purposes, we can drop the
overall constant $N_{c}^{N_{c}}$ that arise when calculating the determinant
of the matrix $W$. Since we are interested in the logarithm of the Haar
measure, we can write:
$\ln\quantity[N_{c}^{N_{c}}\quantity(\ldots)]=N_{c}\ln\quantity[N_{c}]+\ln\quantity[\ldots]$.
The term $N_{c}\ln\quantity[N_{c}]$ can be absorbed in the overall
normalization constant of the generating functional, leaving only the
logarithm of the Haar measure, see Eq. (43). From this definition one can
expect the complexity of the calculation to largely increase for increasing
values of $N_{c}$, leading to polynomials with several interaction terms. For
instance, for $N_{c}=5$, the Haar measure is given by:

Here, $\ell_{\mathrm{F},3}$, $\ell_{\mathrm{F},4}$,
$\overline{\ell}_{\mathrm{F},3}$ and $\overline{\ell}_{\mathrm{F},4}$, are
non-Cartan-Polyakov loops that must be decomposed in terms of Cartan loops
using the previously derived decomposition formulas, in particular, Eqs. (33)
and (34). Written in terms of the Cartan-Polyakov loop set, the Haar measure
for $N_{c}=5$ has 106 different symmetry conserving interaction terms. For
$N_{c}=7$ and $N_{c}=9$ this number increases to 3589 and 164856,
respectively. The reason for this rapid increase in the number of terms can be
traced back to combinatorics: increasing $N_{c}$ increases the number of
elements inside the Cartan set, which in turn, boosts the number of possible
$\mathrm{Z}(N_{c})$ invariant combinations of Cartan-Polyakov loops.

####  IV.1.3 The thermodynamic potential with an arbitrary $N_{c}$

Using the results obtained so far, we are finally able to write a
thermodynamic potential for the PNJL model, in the mean field approximation,
for an arbitrary number of colors:

$\displaystyle\Omega$
$\displaystyle\quantity(N_{c},T,\mu)=U\quantity(N_{c},\expectationvalue{\phi_{i}})+\mathcal{U}_{\mathrm{eff}}\quantity(\mathcal{C}_{N_{c}},T)-2N_{c}\tr\int_{\mathrm{reg}}\frac{\differential[3]{p}}{\quantity(2\pi)^{3}}E$
$\displaystyle-2T\tr\int_{\mathrm{reg}}\frac{\differential[3]{p}}{\quantity(2\pi)^{3}}\ln\quantity[1+\sum_{k=1}^{\frac{N_{c}-1}{2}}c_{k}(-L_{\mathrm{F}})\mathrm{e}^{-k\beta\quantity(E-{\mu})}+\sum_{k=\frac{N_{c}+1}{2}}^{{N_{c}}-1}c_{N_{c}-k}(-L_{\mathrm{F}}^{\dagger})\mathrm{e}^{-k\beta\quantity(E-{\mu})}+\mathrm{e}^{-N_{c}\beta\quantity(E-{\mu})}]$
$\displaystyle-2T\tr\int_{\mathrm{reg}}\frac{\differential[3]{p}}{\quantity(2\pi)^{3}}\ln\quantity[1+\sum_{k=1}^{\frac{N_{c}-1}{2}}c_{k}(-L_{\mathrm{F}}^{\dagger})\mathrm{e}^{-k\beta\quantity(E+{\mu})}+\sum_{k=\frac{N_{c}+1}{2}}^{{N_{c}}-1}c_{N_{c}-k}(-L_{\mathrm{F}})\mathrm{e}^{-k\beta\quantity(E+{\mu})}+\mathrm{e}^{-N_{c}\beta\quantity(E+{\mu})}].$
(73)

As already stated, $\mathcal{C}_{N_{c}}$, is the set of Cartan-Polyakov loops,
the coefficients $c_{k}(-L_{\mathrm{F}})$ and
$c_{k}(-L_{\mathrm{F}}^{\dagger})$ can be calculated with Eq. (22) and the
effective Polyakov loop potential is given by Eq. (60). The mean field
potential, $U\quantity(N_{c},\expectationvalue{\phi_{i}})$, can also be
$N_{c}$ dependent. For instance, in the NJL model, the coupling of the four
scalar-pseudoscalar quark-quark interaction can be considered as being
inversely proportional to $N_{c}$ due to the QCD counting rules [80]. As
already discussed, the values of the mean fields,
$\expectationvalue{\phi_{i}}$, and of the Cartan-Polyakov loops (elements of
the $\mathcal{C}_{N_{c}}$ set) can be obtained by deriving the gap equations
of the model:
$\derivative*{\Omega}{\expectationvalue{\phi_{i}}}=\derivative*{\Omega}{\mathcal{C}_{N_{c}}^{i}}=0$,
where $\mathcal{C}_{N_{c}}^{i}$ is the $i-$th element of the Cartan set.

It is important to note that the fermionic part of the thermodynamic potential
written above, by construction, only depends on Cartan-Polyakov loops.
However, the effective Polyakov loops potential,
$\mathcal{U}_{\mathrm{eff}}\quantity(\mathcal{C}_{N_{c}},T)$, contains the
Haar measure which depends explicitly on non-Cartan-Polyakov loops, in the
fundamental representation, with charges up to $N_{c}-1$ (see Eq. (70)). Of
course, as already discussed, one can make use of the decomposition formula
given in Eq. (20) (or Eq. (25)) in order to write the Haar measure in terms
Cartan-Polyakov loops only. However, when dealing with the mean field
approximation, the values of the Cartan-Polyakov loops are fixed by solving
the gap equations, $\derivative*{\Omega}{\mathcal{C}_{N_{c}}^{i}}=0$. The
derivative of the Haar measure with respect to an element of the Cartan-
Polyakov loop set, $\ell_{\mathrm{F},k}$ (the same can be applied to
derivatives with respect to $\overline{\ell}_{\mathrm{F},k}$), is:

$\displaystyle\derivative{H}{\ell_{\mathrm{F},k}}=H\,\tr\quantity[W\derivative{W}{\ell_{\mathrm{F},k}}].$
(74)

Here, we used Eq. (70) without the $N_{c}^{N_{c}}$ coefficient and
$\derivative*{W}{\ell_{\mathrm{F},k}}$ is a tangent matrix, whose matrix
elements are derivatives of Polyakov loops, in the fundamental representation,
with charges $1,\ldots,N_{c}-1$, with respect to $\ell_{\mathrm{F},k}$ (or
$\overline{\ell}_{\mathrm{F},k}$). Explicitly, it is given by:

When differentiating a Cartan-Polyakov loop, the result is a delta function.
However, the matrix above contains non-Cartan-Polyakov loops with charges in
the interval $(N_{c}+1)/2,\ldots,N_{c}-1$. For these cases one can calculate
these derivatives using:

$\displaystyle\derivative{\ell_{\mathrm{F},a}}{\ell_{\mathrm{F},b}}$
$\displaystyle=\partialderivative{\ell_{\mathrm{F},a}}{\ell_{\mathrm{F},b}}+\sum_{i=\frac{N_{c}+1}{2}}^{a-1}\quantity(\partialderivative{\ell_{\mathrm{F},a}}{\ell_{\mathrm{F},i}})\quantity(\derivative{\ell_{\mathrm{F},i}}{\ell_{\mathrm{F},b}}),$
(75)
$\displaystyle\derivative{\ell_{\mathrm{F},a}}{\overline{\ell}_{\mathrm{F},b}}$
$\displaystyle=\partialderivative{\ell_{\mathrm{F},a}}{\overline{\ell}_{\mathrm{F},b}}+\sum_{i=\frac{N_{c}+1}{2}}^{a-1}\quantity(\partialderivative{\ell_{\mathrm{F},a}}{\ell_{\mathrm{F},i}})\quantity(\derivative{\ell_{\mathrm{F},i}}{\overline{\ell}_{\mathrm{F},b}}),$
(76)
$\displaystyle\derivative{\overline{\ell}_{\mathrm{F},a}}{\ell_{\mathrm{F},b}}$
$\displaystyle=\partialderivative{\overline{\ell}_{\mathrm{F},a}}{\ell_{\mathrm{F},b}}+\sum_{i=\frac{N_{c}+1}{2}}^{a-1}\quantity(\partialderivative{\overline{\ell}_{\mathrm{F},a}}{\overline{\ell}_{\mathrm{F},i}})\quantity(\derivative{\overline{\ell}_{\mathrm{F},i}}{\ell_{\mathrm{F},b}}),$
(77)
$\displaystyle\derivative{\overline{\ell}_{\mathrm{F},a}}{\overline{\ell}_{\mathrm{F},b}}$
$\displaystyle=\partialderivative{\overline{\ell}_{\mathrm{F},a}}{\overline{\ell}_{\mathrm{F},b}}+\sum_{i=\frac{N_{c}+1}{2}}^{a-1}\quantity(\partialderivative{\overline{\ell}_{\mathrm{F},a}}{\overline{\ell}_{\mathrm{F},i}})\quantity(\derivative{\overline{\ell}_{\mathrm{F},i}}{\overline{\ell}_{\mathrm{F},b}}).$
(78)

For more details see Appendix B where all the partial derivatives that are
necessary in the above equations are calculated, for any representation,
$\mathcal{R}$. With this result one can write the gap equations of the PNJL
model for any value of $N_{c}$.

This model can be used, as is, to study the phase diagram of the PNJL model
for increasing values of $N_{c}$. Numerically one can define all the necessary
quantities, impose the gap equations and calculate all thermodynamics
functions of interest for a given value of $N_{c}$. Of course, apart from the
increase in the number of degrees of freedom, increasing $N_{c}$ means that
the largest matrix in $c_{k}(-L_{\mathrm{F}})$ or
$c_{k}(-L_{\mathrm{F}}^{\dagger})$, has dimensions
$(N_{c}-1)/2\times(N_{c}-1)/2$. Thus, naturally, the complexity of the
numerical calculation, using the thermodynamic potential of Eq. (73), gets
larger for larger values of $N_{c}$. Also, Eq. (73) can be a good starting
point to write down a thermodynamic potential for the large $N_{c}$ limit.
Such potential could then be used to study the phase diagram of strongly
interacting matter at large $N_{c}$, and the quarkyonic phase of matter. We
leave both aforementioned studies, the numerical analysis and the
large$-N_{c}$ model, as future research endeavors.

We would like to highlight that the formalism employed here, to define the
PNJL for different values of $N_{c}$, can also be applied to other fermionic
models as long as the thermodynamic potential can be written in a similar way
to the one given in Eq. (44).

###  IV.2 Glue effective potential with an arbitrary $N_{c}$

In this section we apply the concept of Cartan-Polyakov loops to extend the
one-loop gluon effective potential to arbitrary values of $N_{c}$. The
starting point is the usual one-loop contribution to an effective potential of
gluons with a spatially uniform Polyakov loop [81, 82, 57, 83]:

$\displaystyle\Omega_{g}(N_{c},T)=2T\int_{\mathrm{reg}}\frac{\differential[3]{p}}{\quantity(2\pi)^{3}}\quantity{\beta
N_{A}E_{g}+\ln\det\quantity[1-L_{\mathrm{A}}\mathrm{e}^{-\beta E_{g}}]}.$ (79)

Here, $N_{\mathrm{A}}=N_{c}^{2}-1$, is the dimension of the adjoint
representation, $E_{g}=\sqrt{p^{2}+M_{g}^{2}}$ ($M_{g}$ is a phenomenological
gluon effective mass [57, 83]) and, $L_{A}$ is the thermal Wilson line in the
adjoint representation of the $\mathrm{SU}(N_{c})$ group (see Eq. (4)). This
contribution to the effective glue potential drives the system to a state of
spontaneously broken $\mathrm{Z}(N_{c})$ symmetry [44]. In order to study the
transition from a $\mathrm{Z}(N_{c})$ symmetric state to the broken state, one
can consider a confining contribution to the potential in the form of the Haar
measure [44] which, as previously discussed, can be written in terms of the
Cartan-Polyakov loop set, $\mathcal{C}_{N_{c}}$. Hence, the total effective
glue potential, $\Omega(N_{c},T)$, is given by:

$\displaystyle\Omega(N_{c},T)=k(N_{c},T)\ln\quantity[H(\mathcal{C}_{N_{c}})]+\Omega_{g}(N_{c},T).$
(80)

Here, $H(\mathcal{C}_{N_{c}})$ is the Haar measure and $k(N_{c},T)$ is a
phenomenological parameter, similar to the one introduced for the PNJL model
with an arbitrary $N_{c}$ (see the discussion in Section IV.1.2). The
effective glue potential defined above was also derived, for the $N_{c}=3$
case, in Ref. [57] where the background field approximation was applied to the
generating functional of the $\mathrm{SU}(3)$ Yang-Mills theory, in the
presence of a uniform gluon field. This model was also considered in Ref.
[44], in order to study fluctuations of the order parameter in an
$\mathrm{SU}(N_{c})$ effective model.

The strategy to extend this model will be the same as the one used for the
PNJL model. The confining part of the potential, written in Eq. (80), was
previously defined for an arbitrary number of colors (see Eqs. (60) and (70)).
Indeed, in Section IV.1.2 we discussed how to write the Haar measure using
exclusively the elements which constitute the Cartan-Polyakov loop set. The
final step is to evaluate the gluonic determinant, present in Eq. (79), for an
arbitrary number of colors. The gluonic determinant is defined by:

$\displaystyle\det\quantity[1-L_{\mathrm{A}}\mathrm{e}^{-\beta
E_{g}}]=h_{g}^{N_{\mathrm{A}}}\det\quantity[h_{g}^{-1}I-L_{\mathrm{A}}].$ (81)

Like in the case of the fermionic determinant, in the above, we defined the
non-zero quantity, $h_{g}$, as:

$\displaystyle h_{g}=\mathrm{e}^{-\beta E_{g}}.$ (82)

As done for the fermionic determinant, the goal is to calculate the gluonic
determinant for an arbitrary $N_{c}$ using only the corresponding Cartan-
Polyakov loop set. To this end, consider the characteristic polynomial
$p\quantity(L_{\mathrm{A}},h_{g}^{-1})$ (see Eq. (10)). We can write the
gluonic determinant as:

$\displaystyle
h_{g}^{N_{\mathrm{A}}}\det\quantity[h_{g}^{-1}I-L_{\mathrm{A}}]$
$\displaystyle=1+\sum_{k=1}^{N_{c}^{2}-2}c_{k}(L_{\mathrm{A}})h_{g}^{k}+\quantity(-1)^{(N_{c}^{2}-1)}h_{g}^{(N_{c}^{2}-1)}.$
(83)

Where we used the fact that $\det\quantity[L_{\mathrm{A}}]=1$. The
coefficients, $c_{k}(L_{\mathrm{A}})$, which arise in the sum, are calculated
using Eq. (11) with the appropriate choice of representation. Observing the
structure of the matrix $J_{\mathrm{A},k}$, defined in Eq.
(LABEL:JRepk_matrix_def), one can conclude that the $k$ coefficient will
depend on the non-Cartan-Polyakov loops $\ell_{\mathrm{A},k}$, with charges
$1,\ldots,k$. Since we want to write the model in terms of Cartan-Polyakov
loops, we must find a way to write the non-Cartan-Polyakov loop of charges
$k$, in the adjoint representation, $\ell_{\mathrm{A},k}$, in terms of Cartan-
Polyakov loops. In order to accomplish this one can use some results from
group theory.

A known result from representation theory allows for the decomposition of the
product between the fundamental and the conjugate representations of the
$\mathrm{SU}(N_{c})$ group ($\mathrm{F}$ and $\overline{\mathrm{F}}$,
respectively) with its adjoint and the singlet representations [84]:

$\displaystyle\mathrm{F}\otimes\overline{\mathrm{F}}=\mathrm{A}\oplus 1.$ (84)

This result can be proven using Young tableaux or by sandwiching the
fundamental and anti-fundamental thermal Wilson lines between two
$\mathrm{SU}(N_{c})$ generators in order to get an object with adjoint indexes
[39]. Since powers of the thermal Wilson line, $L^{k}$, are elements of the
$\mathrm{SU}(N_{c})$ group, we can write the following identities using
character theory [83]:

$\displaystyle\tr\quantity[L_{\rho\oplus\sigma}^{k}]$
$\displaystyle=\tr\quantity[L_{\rho}^{k}]+\tr\quantity[L_{\sigma}^{k}],$ (85)
$\displaystyle\tr\quantity[L_{\rho\otimes\sigma}^{k}]$
$\displaystyle=\tr\quantity[L_{\rho}^{k}]\tr\quantity[L_{\sigma}^{k}].$ (86)

Where $\rho$ and $\sigma$ are different representations of the
$\mathrm{SU}(N_{c})$ group. Using these properties, alongside the
decomposition made in Eq. (84), one can write:

$\displaystyle\tr[L_{\mathrm{F}\otimes\overline{\mathrm{F}}}^{k}\big{]}$
$\displaystyle=\tr\quantity[L_{\mathrm{A}\oplus 1}^{k}]\Leftrightarrow$
$\displaystyle\tr[L_{\mathrm{F}}^{k}\big{]}\tr[{L_{\mathrm{F}}^{\dagger}}^{k}\big{]}$
$\displaystyle=\tr[L_{\mathrm{A}}^{k}\big{]}+1\Leftrightarrow$ $\displaystyle
N_{\mathrm{F}}^{2}\ell_{\mathrm{F},k}\overline{\ell}_{\mathrm{F},k}$
$\displaystyle=N_{\mathrm{A}}\ell_{\mathrm{A},k}+1.$ (87)

Thus, one can write the Polyakov loop with charge $k$, in the adjoint
representation, in terms of Polyakov loops with charge $k$, in the fundamental
and antifundamental representations. Manipulating the equation above yields
[39, 34]:

$\displaystyle\ell_{\mathrm{A},k}$
$\displaystyle=\frac{1}{1-1/N_{c}^{2}}\ell_{\mathrm{F},k}\overline{\ell}_{\mathrm{F},k}-\frac{1}{N_{c}^{2}-1}.$
(88)

As mentioned for the non-Cartan loops in the fundamental representation, this
relation can also be checked by hand, by considering a particular gauge
choice, e.g., using the already mentioned Polyakov loop gauge [57].

With this result, one can readily write the gluonic determinant as a function
of fundamental and anti-fundamental Polyakov loops of different charges.
Moreover, recursively, one can use the non-Cartan-Polyakov loop decomposition
formula of Eq. (20) for loops with charge $k<N_{c}$ and Eq. (28) for loops
with charge $k\geq N_{c}$, in order to write the gluonic determinant of Eq.
(83) in terms of Cartan-Polyakov loops only. We point out that, just like for
the derivation of the non-Cartan decomposition formulas, or the calculation of
fermionic determinant, one can make use of the identity between the
characteristic polynomial coefficients of unitary matrices, written in Eq.
(18) and shown in Appendix A. Making use of this identity implies that only
$(N_{c}^{2}-1)/2$ different coefficients have to be calculated, with the
remaining ones arising automatically. For $N_{c}=3$, there are a total of
seven coefficients, but only four have to be calculated. Explicitly, for
$N_{c}=3$, one gets:

$\displaystyle c_{1}(L_{\mathrm{A}})=c_{7}(L_{\mathrm{A}})$
$\displaystyle=-8\ell_{\mathrm{A},1}$
$\displaystyle=1-9\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1},$ (89)
$\displaystyle c_{2}(L_{\mathrm{A}})=c_{6}(L_{\mathrm{A}})$
$\displaystyle=32\ell_{\mathrm{A},1}^{2}-4\ell_{\mathrm{A},2}$
$\displaystyle=1-27(\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}-\ell_{\mathrm{F},1}^{3}-\overline{\ell}_{\mathrm{F},1}^{3}),$
(90) $\displaystyle c_{3}(L_{\mathrm{A}})=c_{5}(L_{\mathrm{A}})$
$\displaystyle=-\frac{8}{3}\quantity(32\ell_{\mathrm{A},1}^{3}-12\ell_{\mathrm{A},1}\ell_{\mathrm{A},2}+\ell_{\mathrm{A},3})$
$\displaystyle=-81\ell_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},1}^{2}+27\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}-2,$
(91) $\displaystyle c_{4}(L_{\mathrm{A}})$
$\displaystyle=\frac{8}{3}\quantity(64\ell_{\mathrm{A},1}^{4}-48\ell_{\mathrm{A},2}\ell_{\mathrm{A},1}^{2}+8\ell_{\mathrm{A},1}\ell_{\mathrm{A},3}+3\ell_{\mathrm{A},2}^{2})-2\ell_{\mathrm{A},4}$
$\displaystyle=162\ell_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},1}^{2}+18\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}-54\overline{\ell}_{\mathrm{F},1}^{3}-54\ell_{\mathrm{F},1}^{3}-2.$
(92)

Here, we have written the coefficients first in terms of charged Polyakov
loops in the adjoint representation and then used Eq. (88) alongside the non-
Cartan decomposition formulas, in order to write the coefficients in terms of
Cartan-Polyakov loops. These results agree with the ones present in the
literature, see Ref. [57]. For $N_{c}=5$, one has to calculate twelve
coefficients, the first three are given by:

$\displaystyle c_{1}(L_{\mathrm{A}})=c_{23}(L_{\mathrm{A}})$
$\displaystyle=-24\ell_{\mathrm{A},1}$
$\displaystyle=1-25\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1},$ (93)
$\displaystyle c_{2}(L_{\mathrm{A}})=c_{22}(L_{\mathrm{A}})$
$\displaystyle=288\ell_{\mathrm{A},1}^{2}-12\ell_{\mathrm{A},2}$
$\displaystyle=1+\frac{625}{2}\ell_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},1}^{2}-25\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}-\frac{25}{2}\ell_{\mathrm{F},2}\overline{\ell}_{\mathrm{F},2},$
(94) $\displaystyle c_{3}(L_{\mathrm{A}})=c_{21}(L_{\mathrm{A}})$
$\displaystyle=-8\quantity(288\ell_{\mathrm{A},1}^{3}-36\ell_{\mathrm{A},1}\ell_{\mathrm{A},2}+\ell_{\mathrm{A},3})$
$\displaystyle=1-\frac{15625}{4}\ell_{\mathrm{F},1}^{3}\overline{\ell}_{\mathrm{F},1}^{3}+\frac{3125}{4}\ell_{\mathrm{F},1}^{3}\overline{\ell}_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},2}-\frac{625}{4}\ell_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},1}^{2}+\frac{375}{4}\ell_{\mathrm{F},1}^{2}\overline{\ell}_{\mathrm{F},2}$
$\displaystyle\quad+\frac{3125}{4}\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}\overline{\ell}_{\mathrm{F},1}^{3}-25\ell_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},1}-\frac{625}{4}\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}\overline{\ell}_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},2}+\frac{375}{4}\ell_{\mathrm{F},2}\overline{\ell}_{\mathrm{F},1}^{2}+\frac{375}{4}\overline{\ell}_{\mathrm{F},1}\overline{\ell}_{\mathrm{F},2}^{2}$
$\displaystyle\quad+\frac{3125}{4}\overline{\ell}_{\mathrm{F},1}^{5}-625\overline{\ell}_{\mathrm{F},1}^{3}\overline{\ell}_{\mathrm{F},2}-\frac{125}{4}\ell_{\mathrm{F},2}\overline{\ell}_{\mathrm{F},2}+\frac{3125}{4}\ell_{\mathrm{F},1}^{5}-625\ell_{\mathrm{F},1}^{3}\ell_{\mathrm{F},2}+\frac{375}{4}\ell_{\mathrm{F},1}\ell_{\mathrm{F},2}^{2}.$
(95)

One can immediately observe that the complexity of the coefficients, as
expressions written in terms of Cartan-Polyakov loops, rapidly increases for
larger numbers of colors, as well as, for higher order coefficients for a
fixed $N_{c}$. Indeed, the $c_{12}(L_{\mathrm{A}})$ coefficient, for example,
has 104 different terms when written in terms of the Cartan-Polyakov loop set.
In practice, for numerical calculations, the overall size of the analytical
expressions for a given coefficient at fixed $N_{c}$ is not important. In such
approach, using as variables the Cartan-Polyakov loops, one can build all non-
Cartan-Polyakov loops using the decomposition formulas, make use of Eq. (88)
to directly obtain different charged Polyakov loops in the adjoint
representation and then get the coefficients.

We point out that, since we are dealing with a gluonic theory, the
coefficients only have terms which preserve $\mathrm{Z}(N_{c})$ symmetry.
Naturally, the opposite was found for the fermionic coefficients where the
terms explicitly broke the $\mathrm{Z}(N_{c})$ symmetry.

Finally, one can write the glue effective potential defined in Eq. (79), for
any value of $N_{c}$, using only Cartan-Polyakov loops. It is explicitly given
by:

$\displaystyle\Omega(N_{c},T)$
$\displaystyle=k(N_{c},T)\ln\quantity[H(\mathcal{C}_{N_{c}})]+2\quantity(N_{c}^{2}-1)\int_{\mathrm{reg}}\frac{\differential[3]{p}}{\quantity(2\pi)^{3}}E_{g}$
$\displaystyle+2T\int_{\mathrm{reg}}\frac{\differential[3]{p}}{\quantity(2\pi)^{3}}\ln\quantity[1+\sum_{k=1}^{N_{c}^{2}-2}c_{k}(L_{\mathrm{A}})\mathrm{e}^{-k\beta
E_{g}}+\quantity(-1)^{(N_{c}^{2}-1)}\mathrm{e}^{-(N_{c}^{2}-1)\beta E_{g}}].$
(96)

It is known that the mean field approximation is not suitable to study gluon
effective potentials written in terms of fundamental and anti-fundamental
Polyakov loops [78]. The reason for such can be understood in the basis of the
relation between the adjoint and the fundamental loops, written in Eq. (88).
In the naive mean field approximation, considered in this work, one can see
that in the confined phase, where the expectation values of charged Polyakov
loops in the fundamental representation are zero, the expectation values of
charged adjoint Polyakov loops are negative [78]. As pointed out in Ref. [34],
this implies that charged Polyakov loops in the adjoint representation,
$\ell_{\mathrm{A},k}$, cannot be used as an order parameters for gluon
confinement. Also, the thermodynamic potential obtained in Eq. (96) leads to a
negative pressure, negative entropy density and thus, to thermodynamic
instability [78, 57]. This was discussed in Ref. [57], using the glue model of
Eq. (79) for $N_{c}=3$, by considering the low temperature expansion of the
effective action. In such limit, the last term inside the logarithm,
$\quantity(-1)^{(N_{c}^{2}-1)}\exp\quantity[-(N_{c}^{2}-1)\beta E_{g}]$,
prevails and the incorrect sign for Bose statistics in front of the
exponential is obtained for $N_{c}=3$ [57]. To overcome this feature of the
glue effective potential for $N_{c}=3$, an hybrid model of glueballs in the
confined phase and Polyakov loops in the deconfined phase was suggested, with
the transition between models occurring at the transition temperature. For
further details see Ref. [57]. We point out an interesting feature of the
arbitrary $N_{c}$ effective potential written in Eq. (96): the incorrect Bose
statistics sign is obtained for and odd number of colors while, for an even
number of colors, the sign is the one expected for a system of bosons.

Although one can use the approach used here to extended glue effective
potentials to arbitrary values of $N_{c}$, as discussed above, the simple mean
field approach might not be sufficient and corrections coming from connected
diagrams in Eq. (88), might be necessary [34]. However, at large $N_{c}$,
these corrections become negligible since expectation values factorize [39].
Thus, in this limit, the relation between the adjoint loops of charge $k$ and
the fundamental and anti-fundamental loops of charge $k$ become much simpler:

$\displaystyle\ell_{\mathrm{A},k}$
$\displaystyle\sim\ell_{\mathrm{F},k}\overline{\ell}_{\mathrm{F},k}+O\quantity(\frac{1}{N_{c}}).$
(97)

In fact, at large $N_{c}$, the expectation value of any Polyakov loop can be
written as powers of the expectation value of the fundamental and anti-
fundamental Polyakov loops [39]. Hence, at large $N_{c}$, this model is a
proper tool to study the behavior of the deconfining transition of strongly
interacting matter.

As a final note, we point out that this model can be used to extend the PNJL
model with an arbitrary $N_{c}$ obtained earlier in this work. For instance,
instead of using the Haar measure and the kinetic term proposed in Eq. (59)
for the effective Polyakov loop potential in the PNJL model, one could use the
above effective glue potential. Indeed, one can show that, for $N_{c}=3$,
under certain assumptions, it is possible to recover from Eq. (96) the
$-\overline{\Phi}\Phi$ used in the usual $N_{c}=3$ PNJL model, which is
essential in order to get a first-order phase transition [57]. Not only that,
using this effective glue potential, for the gluonic sector of PNJL model,
would remove the necessity of defining so many $N_{c}-$dependent
phenomenological parameters present in the effective Polyakov loop potential
(see Eq. (60)). In such case, the only parameter coming from the gluonic
sector would be the one in front of the Haar measure [57]. Of course, the
problem of the negative contribution to the pressure in the confined phase,
for $\ell_{\mathrm{F},k}=0$, would have to be addressed [57, 44]. This issue
is partially resolved in Ref. [44] by incorporating, in the thermodynamic
potential, a contribution coming from the ghost fields. Such contribution
reinforces the confining part of the potential, opposing the gluonic
(deconfining) part and it might be able to overturn the negative tendency of
the gluonic pressure in the confined phase [44]. In the same work, it is
discussed a phenomenological solution to this problem for $N_{c}=3$: subtract,
from the glue effective potential, its contribution at
$\ell_{\mathrm{F},1}=\overline{\ell}_{\mathrm{F},1}=0$. Employing this recipe
does not change the results for the Polyakov loops, but requires one to
correct the number of gluonic degrees of freedom at large temperatures.
Naturally, for other values of $N_{c}$, this phenomenological approach would
also work, albeit one would have to subtract the potential at
$\ell_{\mathrm{F},k}=\overline{\ell}_{\mathrm{F},k}=0$, with
$k=1,\ldots,(N_{c}-1)/2$, i.e., considering all elements of the Cartan-
Polyakov loop set to vanish. To explore additional works that address and
offer solutions to this problem, beyond the ones already mentioned, see Refs.
[85, 86, 87, 88].

##  V Conclusions

In this work we defined Cartan-Polyakov loops, a special subset of Polyakov
loops in the fundamental representation of the $\mathrm{SU}(N_{c})$ group
which, under $\mathrm{Z}(N_{c})$ transformations, behave like fields with
charges $k=1,\ldots,(N_{c}-1)/2$. This set can be used as independent degrees
of freedom to parametrize effective models of quarks and gluons with any
number of colors. When building such phenomenological models, non-Cartan-
Polyakov loops, Polyakov loops with charges higher than $(N_{c}-1)/2$,
naturally arise. In order to write these models as functions of the loops
contained in the Cartan set, we derived a decomposition formula for non-
Cartan-Polyakov loops. Such formula was derived using properties of the
characteristic polynomial of unitary matrices and the Cayley-Hamilton theorem.
Fixing some values of $N_{c}$, we consistently decomposed some non-Cartan
loops until they were only functions of Cartan-Polyakov loops. Finally, we
showed how to build two distinct effective models for any $N_{c}$, including a
version of the PNJL model and an effective glue potential. These models were
built using the naive mean field approximation commonly used when dealing with
the PNJL model: the Cartan-Polyakov loops were directly substituted by their
respective expectation values and the expectation values of products between
them factorized. We highlight that the potentials presented in this work, were
derived using a particular gauge choice and/or required some phenomenological
parameters (see Eq. (73) for the thermodynamic potential of the PNJL model and
Eq. (96) for the glue effective potential, both written for an arbitrary
$N_{c}$). For instance, the number of colors and temperature dependent
prefactors $B\quantity(N_{c},T)$ and $k(N_{c},T)$, in Eqs. (60) and (80)
respectively, were included because of phenomenological considerations. They
were not derived from pure gauge theory.

The formalism developed in this work, allows for a consistent study on the
effect of increasing the number of colors in effective models of quarks and
gluons. Numerically, one can build tools that, for a given $N_{c}$,
automatically builds the Cartan-Polyakov loop set and calculate all the
necessary non-Cartan-Polyakov loops from it. In the case of the PNJL model,
for example, one can build numerical tools in order to solve the model, within
the mean field approximation, for a fixed value of $N_{c}$. Having the
solution of such model, one can analyze the phase diagram of the model for
that particular value of $N_{c}$, especially the behavior of the critical
endpoint and the interplay between the chiral transition and the deconfinement
transition. One can then repeat the calculation to increasing values of
$N_{c}$ and study the behavior of the critical endpoint for large values of
$N_{c}$, as well as study the conjectured quarkyonic phase of matter at high
chemical potentials and moderate temperatures. Indeed, in Ref. [29], the fate
of the critical endpoint at large $N_{c}$ was studied using a Polyakov loop
quark-meson model by varying the number of colors. Within this work it was
also confirmed the existence of the quarkyonic phase, in which chiral symmetry
is restored but quarks are still confined [29]. It would be very interesting
to see if similar features were obtained within the model proposed in this
work. In the case of the glue effective model defined in this work, it can be
used to study the fluctuations of the confined-deconfined order parameters,
given by different susceptibilities of the Polyakov loops, and the curvature
masses associated with the Cartan-Polyakov loop set, for any value of $N_{c}$.
These observables can be extremely useful tools to study the deconfinement
phase transition in the Yang-Mills theory [89], providing useful information
about the link between the structure of the vacuum and gluon properties [44].
For example, in Ref. [89], the susceptibilities of the modulus, the real part
and the imaginary part of the Polyakov loop, in the fundamental
representation, were calculated in $\mathrm{SU}(3)$ lattice gauge theory (for
studies including quark fields see Refs. [73, 90, 91]). In that work the
authors were able to show that different ratios of these susceptibilities are
independent of the renormalization of the Polyakov and that these ratios
depend weakly on the system size. Additionally, from the point of view of the
phase transition, in the same work, it was shown that these ratios are
discontinuous at the phase transition point and that above and below this
point they are almost temperature independent. For other studies involving the
susceptibilities of the Polyakov loop, see Refs. [73, 92, 93]. A similar study
carried out for larger values of $N_{c}$, using $\mathrm{SU}(N_{c})$ lattice
gauge theory, might be extremely useful in order to verify the relevance of
the Cartan-Polyakov loop set. In such study one can evaluate different
observables in the form of different ratios of different elements of the
Cartan-Polyakov loop set. For instance, for $N_{c}=5$, the Cartan-Polyakov set
is
$\mathcal{C}_{5}=\quantity{\ell_{\mathrm{F},1},\ell_{\mathrm{F},2},\overline{\ell}_{\mathrm{F},1},\overline{\ell}_{\mathrm{F},2}}$
and the ratio of the modulus, and real and imaginary parts of distinct
elements of the set, as a function of temperature, can be analyzed.

The Cartan-Polyakov loop set can also be relevant to study the so-called
Casimir scaling hypothesis. For Polyakov loops, this hypothesis states that
there is a scaling relation between the Polyakov loops in different
representations of the $\mathrm{SU}(N_{c})$ group [40, 41, 83]. The Casimir
hypothesis is not only applicable to Polyakov loops but also to other
observables [94]. For values of $N_{c}>3$, it can be insightful to understand
if the hypothesis holds for all elements of the Cartan-Polyakov loop set,
besides the $\ell_{\mathrm{F},1}$ and $\overline{\ell}_{\mathrm{F},1}$. For
some works where the Casimir scaling hypothesis is discussed and studied, see
Refs. [95, 40, 41, 83, 96, 96, 94].

As future work, we plan to use the formalism developed in this paper to build
numerical tools which are prepared to solve the PNJL model for any given
$N_{c}$, with the goal of studying the behavior of the critical endpoint with
increasing $N_{c}$. As pointed out, this can be accomplished recursively
without further approximations, except the ones introduced in this work. At
the same time, we will build an effective PNJL model at large $N_{c}$. In this
limit, we plan to study the phase diagram of the model and the properties of
the expected quarkyonic phase at large baryon densities. Having a numerical
calculation together with a large $N_{c}$ approximation will allow us to
compare both approaches and better understand the behavior of strongly
interacting matter at large $N_{c}$.

## Acknowledgments

The authors would like to thank João Moreira and Orlando Oliveira for useful
discussions. The authors also thank Tomáš Brauner for elucidating how to
express the $\mathrm{SU}(3)$ Haar measure in terms of the character of the
fundamental representation, in an online document entitled “Haar measure on
the unitary groups” [79]. This work was supported by project CERN/FIS-
PAR/0040/2019.

##  Appendix A The characteristic polynomial

The characteristic polynomial of a square $N\times N$ matrix $A$ is a
polynomial, $p\quantity(A,\lambda)$, in the variable $\lambda$, whose roots
are the eigenvalues of the matrix $A$. This polynomial is defined by [97, 98,
99]:

$\displaystyle p\quantity(A,\lambda)\equiv\det\quantity[A-\lambda I],$ (98)

where $I$ is the identity matrix. It can be shown that the characteristic
polynomial can be written as:

$\displaystyle
p\quantity(A,\lambda)=\quantity(-1)^{N}\quantity[\lambda^{N}+c_{1}(A)\lambda^{N-1}+c_{2}(A)\lambda^{N-2}+\dots+c_{N-1}(A)\lambda+c_{N}(A)].$
(99)

By considering the special case of $\lambda=0$, the $c_{N}$ coefficient is
readily obtained. It is given by:

$\displaystyle c_{N}(A)=\quantity(-1)^{N}\det\quantity[A].$ (100)

The coefficients $c_{k}(A)$ can be calculated using [97, 98, 99]: Here,
$t_{k}=\tr\quantity[A^{k}]$. One can also formally define a zeroth order
coefficient, $c_{0}=1$.

###  A.1  Relation between coefficients for $\mathrm{U}(N)$ matrices

Consider an unitary matrix $M$ with dimension $N\times N$. Using Eq. (99), its
characteristic polynomial, $p\quantity(M,\lambda)$, is:

$\displaystyle p\quantity(M,\lambda)$ $\displaystyle=\det\quantity[M-\lambda
I]$
$\displaystyle=\quantity(-1)^{N}\quantity(\lambda^{N}+\sum_{i=1}^{N-1}c_{i}(M)\lambda^{N-i}+\quantity(-1)^{N}\det\quantity[M]).$
(101)

Using this definition and the properties of the unitary matrix $M$, one can
write the polynomial $p\quantity(M,\lambda)$ as follows:

$\displaystyle p\quantity(M,\lambda)$ $\displaystyle=\det\quantity[M-\lambda
MM^{\dagger}]$ $\displaystyle=\det\quantity[-\lambda
M\quantity(M^{\dagger}-\overline{\lambda}I)]$
$\displaystyle=\det\quantity[-\lambda
M]\det\quantity[M^{\dagger}-\overline{\lambda}I]$
$\displaystyle=\quantity(-1)^{N}\lambda^{N}\det\quantity[M]\det\quantity[M^{\dagger}-\overline{\lambda}I].$
(102)

Here, we assumed that $\lambda$ is non-zero and defined
$\overline{\lambda}=\lambda^{-1}$. The determinant on the right hand side of
this equation, $\det\quantity[M^{\dagger}-\overline{\lambda}I]$, can also be
written in terms of a characteristic polynomial,
$p\quantity(M^{\dagger},\overline{\lambda})$. It yields:

$\displaystyle p\quantity(M,\lambda)$
$\displaystyle=\quantity(-1)^{N}\lambda^{N}\det\quantity[M]\,p\quantity(M^{\dagger},\overline{\lambda})$
$\displaystyle=\quantity(-1)^{N}\lambda^{N}\det\quantity[M]\quantity(\quantity(-1)^{N}\overline{\lambda}^{N}+\quantity(-1)^{N}\sum_{i=1}^{N-1}c_{i}(M^{\dagger})\overline{\lambda}^{N-i}+\det\quantity[M^{\dagger}])$
$\displaystyle=\quantity(-1)^{N}\quantity(\quantity(-1)^{N}\det\quantity[M]+\quantity(-1)^{N}\det\quantity[M]\sum_{i=1}^{N-1}c_{i}(M^{\dagger})\lambda^{i}+\lambda^{N}\det\quantity[MM^{\dagger}])$
$\displaystyle=\quantity(-1)^{N}\quantity(\lambda^{N}+\sum_{i=1}^{N-1}\quantity(-1)^{N}\det\quantity[M]c_{N-i}(M^{\dagger})\lambda^{N-i}+\quantity(-1)^{N}\det\quantity[M]).$
(103)

Comparing this polynomial, with the polynomial in Eq. (101), one can establish
the following identities:

$\displaystyle c_{1}(M)$
$\displaystyle=\quantity(-1)^{N}\det\quantity[M]c_{N-1}(M^{\dagger}),$
$\displaystyle c_{2}(M)$
$\displaystyle=\quantity(-1)^{N}\det\quantity[M]c_{N-2}(M^{\dagger}),$
$\displaystyle\;\;\vdots$ (104) $\displaystyle c_{N-1}(M)$
$\displaystyle=\quantity(-1)^{N}\det\quantity[M]c_{1}(M^{\dagger}).$

For a general coefficient with index $j$, one can write:

$\displaystyle c_{j}(M)$
$\displaystyle=\quantity(-1)^{N}\det\quantity[M]c_{N-j}(M^{\dagger}),$ (105)

with $j=0,1,2,3,\dots,N-1,N$. We extended the relation for the trivial $0-$th
and $N-$th order coefficients, $c_{0}(M)$ and $c_{N}(M)$. This extension can
be checked explicitly (by definition, $c_{0}=1$):

$\displaystyle c_{0}(M)$
$\displaystyle=\quantity(-1)^{N}\det\quantity[M]c_{N}(M^{\dagger})=\quantity(-1)^{2N}\det\quantity[MM^{\dagger}]=1.$

In the above we made use of Eq. (100). We highlight that $N$ is the dimension
of the unitary matrix $M$ and $c_{j}(M)$ is calculated from Eq.
(LABEL:c_k_definition).

##  Appendix B Differentiating non-Cartan-Polyakov loops

For a given $N_{c}$, one can write a non-Cartan-Polyakov loop of charge $k$ in
the $\mathcal{R}$ representation of the $\mathrm{SU}(N_{c})$ group, as a
function of the Cartan-Polyakov loop set for that specific $N_{c}$, using the
decomposition formulas introduced in Eqs.(20), (25) and (28). For more
details, see Section III. Hence, one can write a formula for the total
derivative of a non-Cartan-Polyakov loop with respect to some element of the
Cartan set. This result is very important in the evaluation of the gap
equations of the PNJL model, arising when one is calculating derivatives of
the Haar measure with respect to Cartan-Polyakov loops.

Consider a non-Cartan-Polyakov loop of charge $a$, in some representation of
the $\mathrm{SU}(N_{c})$ group, $\ell_{\mathcal{R},a}$. As discussed in this
work, one can always build these quantities recursively using Cartan loops. In
practice, this non-Cartan loop can be considered as function of the elements
in the Cartan set and as a function of all non-Cartan loops of smaller charge,
which, in turn, are functions of the Cartan set and so on. Thus, we can write:

$\displaystyle\ell_{\mathcal{R},a}=\ell_{\mathcal{R},a}\quantity(\mathcal{C}_{N_{c}},\ell_{\mathcal{R},a-1},\ell_{\mathcal{R},a-2},\ldots).$
(106)

Thus, the total derivative of this non-Cartan loop, $\ell_{\mathcal{R},a}$,
with respect to a Cartan-Polyakov loop $\ell_{\mathcal{R},b}$, is given by the
following expression:

$\displaystyle\derivative{\ell_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}$
$\displaystyle=\partialderivative{\ell_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}+\sum_{i=\frac{N_{c}+1}{2}}^{a-1}\quantity(\partialderivative{\ell_{\mathcal{R},a}}{\ell_{\mathcal{R},i}})\quantity(\derivative{\ell_{\mathcal{R},i}}{\ell_{\mathcal{R},b}}).$
(107)

This formula also works recursively: in order to get
$\derivative*{\ell_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}$ one needs to first
calculate all the total derivatives of non-Cartan loops with smaller charges,
$\derivative*{\ell_{\mathcal{R},i}}{\ell_{\mathcal{R},b}}$, with
$i\in\quantity{\nicefrac{{(N_{c}+1)}}{{2}},\ldots,a-1}$. Likewise, one can
write the total derivative of a non-Cartan loop with respect to some adjoint
Cartan loop,
$\derivative*{\ell_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}$ or
$\derivative*{\overline{\ell}_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}$ and
$\derivative*{\overline{\ell}_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}$.
They are given by:

$\displaystyle\derivative{\ell_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}$
$\displaystyle=\partialderivative{\ell_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}+\sum_{i=\frac{N_{c}+1}{2}}^{a-1}\quantity(\partialderivative{\ell_{\mathcal{R},a}}{\ell_{\mathcal{R},i}})\quantity(\derivative{\ell_{\mathcal{R},i}}{\overline{\ell}_{\mathcal{R},b}}),$
(108)
$\displaystyle\derivative{\overline{\ell}_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}$
$\displaystyle=\partialderivative{\overline{\ell}_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}+\sum_{i=\frac{N_{c}+1}{2}}^{a-1}\quantity(\partialderivative{\overline{\ell}_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},i}})\quantity(\derivative{\overline{\ell}_{\mathcal{R},i}}{\ell_{\mathcal{R},b}}),$
(109)
$\displaystyle\derivative{\overline{\ell}_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}$
$\displaystyle=\partialderivative{\overline{\ell}_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}+\sum_{i=\frac{N_{c}+1}{2}}^{a-1}\quantity(\partialderivative{\overline{\ell}_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},i}})\quantity(\derivative{\overline{\ell}_{\mathcal{R},i}}{\overline{\ell}_{\mathcal{R},b}}).$
(110)

The different partial derivatives, necessary to evaluate the total
derivatives, can be obtained by partially differentiating the decomposition
formulas derived in this work. If the charge $a$ is smaller than
$N_{\mathcal{R}}$, the partial derivative can be calculated using the
decomposition formula of Eq. (20) (a similar calculation can be done using Eq.
(25) as a starting point) while, for charges $a\geq N_{\mathcal{R}}$, one can
calculate the partial derivative using Eq. (28). Lets first deal with the case
$a<N_{\mathcal{R}}$. Calculating the partial derivative of Eq. (20), with
respect to the Cartan-Polyakov loop $\ell_{\mathcal{R},b}$, yields:

$\displaystyle\left.\partialderivative{\ell_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}\right|_{a<N_{\mathcal{R}}}$
$\displaystyle=(-1)^{a+1}\frac{a}{N_{\mathcal{R}}\quantity(N_{\mathcal{R}}-a)!}\partialderivative{\ell_{\mathcal{R},b}}\quantity(\det\quantity[\overline{J}_{\mathcal{R},N_{\mathcal{R}}-a}])-\frac{a}{a!}\sum_{i=1}^{a-1}\quantity(-1)^{i-a}\partialderivative{\ell_{\mathcal{R},b}}\quantity(\ell_{\mathcal{R},i}\det\quantity[J_{\mathcal{R},a}^{(i,1)}])$
$\displaystyle=-\frac{a}{a!}\quantity(\quantity(-1)^{b-a}\det\quantity[J_{\mathcal{R},a}^{(b,1)}]+\sum_{i=1}^{a-1}\quantity(-1)^{i-a}\ell_{\mathcal{R},i}\det\quantity[J_{\mathcal{R},a}^{(i,1)}]\tr\quantity[\quantity[J_{\mathcal{R},a}^{(i,1)}]^{-1}\partialderivative{J_{\mathcal{R},a}^{(i,1)}}{\ell_{\mathcal{R},b}}]).$
(111)

Here, we used Jacobi’s formula in order to express the derivative of the
determinant, $\derivative{t}\det[A(t)]=$ $\det[A(t)]\times$
$\tr[A(t)^{-1}\derivative{A(t)}{t}]$. We also used the fact that
$\partialderivative{\ell_{\mathcal{R},b}}\quantity(\det\quantity[\overline{J}_{\mathcal{R},N_{\mathcal{R}}-a}])=0$,
since the matrix $\overline{J}_{\mathcal{R},N_{\mathcal{R}}-a}$ only depends
on adjoint Cartan-Polyakov loops. We recall that in the notation used
throughout this work, the subscript $(i,1)$ in the matrix
$J_{\mathcal{R},a}^{(i,1)}$ is equivalent to the matrix $J_{\mathcal{R},a}$
with the $i-$th line and first column removed.

The matrix $\partialderivative{J_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}$ can
be found very easily by taking into account its definition given in Eq.
(LABEL:JRepk_matrix_def). It is explicitly given by:

If the non-Cartan loop has charge $a\geq N_{\mathcal{R}}$, one must calculate
the partial derivative of Eq. (20). It yields:

$\displaystyle\left.\partialderivative{\ell_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}\right|_{a\geq
N_{\mathcal{R}}}$
$\displaystyle=\quantity(-1)^{N_{\mathcal{R}}+1}\partialderivative{\ell_{\mathcal{R},b}}\ell_{\mathcal{R},a-N_{\mathcal{R}}}-\sum_{i=1}^{{N_{\mathcal{R}}}-1}\partialderivative{\ell_{\mathcal{R},b}}\quantity(c_{i}(L_{\mathcal{R}})\ell_{\mathcal{R},a-i})$
$\displaystyle=\quantity(-1)^{N_{\mathcal{R}}+1}\delta_{b,a-N_{\mathcal{R}}}-\sum_{i=1}^{{N_{\mathcal{R}}}-1}\delta_{b,a-i}c_{i}(L_{\mathcal{R}})-\sum_{i=1}^{{N_{\mathcal{R}}}-1}c_{i}(L_{\mathcal{R}})\ell_{\mathcal{R},a-i}\tr\quantity[\quantity[J_{\mathcal{R},a}]^{-1}\partialderivative{J_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}].$
(113)

Likewise, one can calculate the other partial derivatives. Summarizing, one
gets:

$\displaystyle\left.\partialderivative{\overline{\ell}_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}\right|_{a<N_{\mathcal{R}}}$
$\displaystyle=-\frac{a}{a!}\quantity(\quantity(-1)^{b-a}\det\quantity[\overline{J}_{\mathcal{R},a}^{(b,1)}]+\sum_{i=1}^{a-1}\quantity(-1)^{i-a}\overline{\ell}_{\mathcal{R},i}\det\quantity[\overline{J}_{\mathcal{R},a}^{(i,1)}]\tr\quantity[\quantity[\overline{J}_{\mathcal{R},a}^{(i,1)}]^{-1}\partialderivative{\overline{J}_{\mathcal{R},a}^{(i,1)}}{\overline{\ell}_{\mathcal{R},b}}]),$
(114)
$\displaystyle\left.\partialderivative{\overline{\ell}_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}\right|_{a\geq
N_{\mathcal{R}}}$
$\displaystyle=\quantity(-1)^{N_{\mathcal{R}}+1}\delta_{b,a-N_{\mathcal{R}}}-\sum_{i=1}^{{N_{\mathcal{R}}}-1}\delta_{b,a-i}c_{i}(L_{\mathcal{R}}^{\dagger})-\sum_{i=1}^{{N_{\mathcal{R}}}-1}c_{i}(L_{\mathcal{R}}^{\dagger})\overline{\ell}_{\mathcal{R},a-i}\tr\quantity[\quantity[\overline{J}_{\mathcal{R},a}]^{-1}\partialderivative{\overline{J}_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}],$
(115)
$\displaystyle\left.\partialderivative{\overline{\ell}_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}\right|_{a<N_{\mathcal{R}}}$
$\displaystyle=(-1)^{a+1}\frac{a}{N_{\mathcal{R}}\quantity(N_{\mathcal{R}}-a)!}\det\quantity[J_{\mathcal{R},N_{\mathcal{R}}-a}]\tr\quantity[\quantity[J_{\mathcal{R},N_{\mathcal{R}}-a}]^{-1}\partialderivative{J_{\mathcal{R},N_{\mathcal{R}}-a}}{\ell_{\mathcal{R},b}}],$
(116)
$\displaystyle\left.\partialderivative{\overline{\ell}_{\mathcal{R},a}}{\ell_{\mathcal{R},b}}\right|_{a\geq
N_{\mathcal{R}}}$ $\displaystyle=0,$ (117)
$\displaystyle\left.\partialderivative{\ell_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}\right|_{a<N_{\mathcal{R}}}$
$\displaystyle=(-1)^{a+1}\frac{a}{N_{\mathcal{R}}\quantity(N_{\mathcal{R}}-a)!}\det\quantity[\overline{J}_{\mathcal{R},N_{\mathcal{R}}-a}]\tr\quantity[\quantity[\overline{J}_{\mathcal{R},N_{\mathcal{R}}-a}]^{-1}\partialderivative{\overline{J}_{\mathcal{R},N_{\mathcal{R}}-a}}{\overline{\ell}_{\mathcal{R},b}}],$
(118)
$\displaystyle\left.\partialderivative{\ell_{\mathcal{R},a}}{\overline{\ell}_{\mathcal{R},b}}\right|_{a\geq
N_{\mathcal{R}}}$ $\displaystyle=0.$ (119)

With these tools one is able to evaluate the total derivatives of quantities
which depend on non-Cartan-Polyakov loops, as is the case with the Haar
measure (see Eq. (70)). So, one can readily evaluate the derivatives of the
Haar measure with respect to the elements of the Cartan-Polyakov loop set and
obtain the gap equations of the PNJL model for any value of $N_{c}$.

## References

  * Skands [2013] P. Skands, Introduction to QCD, in  _Searching for New Physics at Small and Large Scales_ (2013) pp. 341–420, arXiv:1207.2389 [hep-ph] .
  * Halasz  _et al._ [1998] A. M. Halasz, A. D. Jackson, R. E. Shrock, M. A. Stephanov, and J. J. M. Verbaarschot, On the phase diagram of QCD, Phys. Rev. D 58, 096007 (1998), arXiv:hep-ph/9804290 .
  * Hansen  _et al._ [2020] H. Hansen, R. Stiele, and P. Costa, Quark and Polyakov-loop correlations in effective models at zero and nonvanishing density, Phys. Rev. D 101, 094001 (2020), arXiv:1904.08965 [hep-ph] .
  * Schmidt and Sharma [2017] C. Schmidt and S. Sharma, The phase structure of QCD, J. Phys. G 44, 104002 (2017), arXiv:1701.04707 [hep-lat] .
  * ’t Hooft [1974] G. ’t Hooft, A Planar Diagram Theory for Strong Interactions, Nucl. Phys. B 72, 461 (1974).
  * Witten [1979] E. Witten, Baryons in the 1/n Expansion, Nucl. Phys. B 160, 57 (1979).
  * Jenkins [1998] E. E. Jenkins, Large N(c) baryons, Ann. Rev. Nucl. Part. Sci. 48, 81 (1998), arXiv:hep-ph/9803349 .
  * Bonanno and Giacosa [2011] L. Bonanno and F. Giacosa, Does nuclear matter bind at large $N_{c}$?, Nucl. Phys. A 859, 49 (2011), arXiv:1102.3367 [hep-ph] .
  * Lucini and Panero [2014] B. Lucini and M. Panero, Introductory lectures to large-$\scriptsize{N}$ QCD phenomenology and lattice results, Prog. Part. Nucl. Phys. 75, 1 (2014), arXiv:1309.3638 [hep-th] .
  * Richardson  _et al._ [2021] T. R. Richardson, M. R. Schindler, S. Pastore, and R. P. Springer, Large-$N_{c}$ analysis of two-nucleon neutrinoless double-$\beta$ decay and charge-independence-breaking contact terms, Phys. Rev. C 103, 055501 (2021), arXiv:2102.02184 [nucl-th] .
  * Lucha  _et al._ [2021] W. Lucha, D. Melikhov, and H. Sazdjian, Tetraquarks in large-Nc QCD, Prog. Part. Nucl. Phys. 120, 103867 (2021), arXiv:2102.02542 [hep-ph] .
  * Flores-Mendieta  _et al._ [2021] R. Flores-Mendieta, C. I. Garcia, and J. Hernandez, Baryon axial vector current in large-$N_{c}$ chiral perturbation theory: Complete analysis for $N_{c}=3$, Phys. Rev. D 103, 094032 (2021), arXiv:2102.06100 [hep-ph] .
  * Hernández and Romero-López [2021] P. Hernández and F. Romero-López, The large $N_{c}$ limit of QCD on the lattice, Eur. Phys. J. A 57, 52 (2021), arXiv:2012.03331 [hep-lat] .
  * Karthik and Narayanan [2022] N. Karthik and R. Narayanan, Parton physics of the large-Nc mesons, Phys. Rev. D 106, 014503 (2022), arXiv:2205.02252 [hep-lat] .
  * Teper [1997] M. Teper, SU(N(c)) gauge theories for all N(c) in three-dimensions and four-dimensions, Phys. Lett. B 397, 223 (1997), arXiv:hep-lat/9701003 .
  * Teper [1999] M. J. Teper, SU(N) gauge theories in (2+1)-dimensions, Phys. Rev. D 59, 014512 (1999), arXiv:hep-lat/9804008 .
  * Di Vecchia [1999] P. Di Vecchia, Large N gauge theories and AdS / CFT correspondence, in  _ICTP Trieste Spring Workshop on Superstrings and Related Matters_ (1999) pp. 1–70, arXiv:hep-th/9908148 .
  * DeGrand and Liu [2016] T. DeGrand and Y. Liu, Lattice study of large $N_{c}$ QCD, Phys. Rev. D 94, 034506 (2016), [Erratum: Phys.Rev.D 95, 019902 (2017)], arXiv:1606.01277 [hep-lat] .
  * García Pérez [2020] M. García Pérez, Prospects for large N gauge theories on the lattice, PoS LATTICE2019, 276 (2020), arXiv:2001.10859 [hep-lat] .
  * Polchinski [1992] J. Polchinski, High temperature limit of the confining phase, Phys. Rev. Lett. 68, 1267 (1992), arXiv:hep-th/9109007 .
  * Makeenko [2010] Y. Makeenko, A Brief Introduction to Wilson Loops and Large N, Phys. Atom. Nucl. 73, 878 (2010), arXiv:0906.4487 [hep-th] .
  * Maldacena [1998] J. M. Maldacena, The Large N limit of superconformal field theories and supergravity, Adv. Theor. Math. Phys. 2, 231 (1998), arXiv:hep-th/9711200 .
  * Aharony  _et al._ [2000] O. Aharony, S. S. Gubser, J. M. Maldacena, H. Ooguri, and Y. Oz, Large N field theories, string theory and gravity, Phys. Rept. 323, 183 (2000), arXiv:hep-th/9905111 .
  * Yang and Yuan [2022] Y. Yang and P.-H. Yuan, QCD phase diagram by holography, Phys. Lett. B 832, 137212 (2022), arXiv:2011.11941 [hep-th] .
  * McLerran and Pisarski [2007] L. McLerran and R. D. Pisarski, Phases of cold, dense quarks at large N(c), Nucl. Phys. A 796, 83 (2007), arXiv:0706.2191 [hep-ph] .
  * McLerran  _et al._ [2009] L. McLerran, K. Redlich, and C. Sasaki, Quarkyonic Matter and Chiral Symmetry Breaking, Nucl. Phys. A 824, 86 (2009), arXiv:0812.3585 [hep-ph] .
  * Torrieri and Mishustin [2010] G. Torrieri and I. Mishustin, The nuclear liquid-gas phase transition at large $N_{c}$ in the Van der Waals approximation, Phys. Rev. C 82, 055202 (2010), arXiv:1006.2471 [nucl-th] .
  * Buisseret and Lacroix [2012] F. Buisseret and G. Lacroix, A large-$N_{c}$ PNJL model with explicit Z${}_{N_{c}}$ symmetry, Phys. Rev. D 85, 016009 (2012), arXiv:1107.4672 [hep-ph] .
  * Kovács  _et al._ [2022] P. Kovács, G. Kovács, and F. Giacosa, Fate of the critical endpoint at large Nc, Phys. Rev. D 106, 116016 (2022), arXiv:2209.09568 [hep-ph] .
  * Giacosa and Pagliara [2017] F. Giacosa and G. Pagliara, Neutron stars in the large-$N_{c}$ limit, Nucl. Phys. A 968, 366 (2017), arXiv:1707.02644 [nucl-th] .
  * Margueron  _et al._ [2021] J. Margueron, H. Hansen, P. Proust, and G. Chanfray, Quarkyonic stars with isospin-flavor asymmetry, Phys. Rev. C 104, 055803 (2021), arXiv:2103.10209 [nucl-th] .
  * Pisarski [2002] R. D. Pisarski, Notes on the deconfining phase transition, in  _Cargese Summer School on QCD Perspectives on Hot and Dense Matter_ (2002) pp. 353–384, arXiv:hep-ph/0203271 .
  * Fukushima [2004] K. Fukushima, Chiral effective model with the Polyakov loop, Phys. Lett. B591, 277 (2004), arXiv:hep-ph/0310121 [hep-ph] .
  * Fukushima and Skokov [2017] K. Fukushima and V. Skokov, Polyakov loop modeling for hot QCD, Prog. Part. Nucl. Phys. 96, 154 (2017), arXiv:1705.00718 [hep-ph] .
  * McLerran and Svetitsky [1981] L. D. McLerran and B. Svetitsky, Quark Liberation at High Temperature: A Monte Carlo Study of SU(2) Gauge Theory, Phys. Rev. D 24, 450 (1981).
  * [36] C. Ratti, Thermodynamics of the quark-gluon plasma, in  _International School on Quark-Gluon Plasma and Heavy Ion Collisions : past, present, future_.
  * Zee [2016] A. Zee,  _Group Theory in a Nutshell for Physicists_ (Princeton University Press, USA, 2016).
  * Hidaka and Pisarski [2010] Y. Hidaka and R. D. Pisarski, Small shear viscosity in the semi quark gluon plasma, Phys. Rev. D 81, 076002 (2010), arXiv:0912.0940 [hep-ph] .
  * Dumitru  _et al._ [2004] A. Dumitru, Y. Hatta, J. Lenaghan, K. Orginos, and R. D. Pisarski, Deconfining phase transition as a matrix model of renormalized Polyakov loops, Phys. Rev. D 70, 034511 (2004), arXiv:hep-th/0311223 .
  * Gupta  _et al._ [2008] S. Gupta, K. Huebner, and O. Kaczmarek, Renormalized Polyakov loops in many representations, Phys. Rev. D 77, 034503 (2008), arXiv:0711.2251 [hep-lat] .
  * Abuki and Fukushima [2009] H. Abuki and K. Fukushima, Gauge dynamics in the PNJL model: Color neutrality and Casimir scaling, Phys. Lett. B 676, 57 (2009), arXiv:0901.4821 [hep-ph] .
  * Ayyar  _et al._ [2018] V. Ayyar, T. DeGrand, D. C. Hackett, W. I. Jay, E. T. Neil, Y. Shamir, and B. Svetitsky, Finite-temperature phase structure of SU(4) gauge theory with multiple fermion representations, Phys. Rev. D 97, 114502 (2018), arXiv:1802.09644 [hep-lat] .
  * Megias  _et al._ [2006] E. Megias, E. Ruiz Arriola, and L. L. Salcedo, Polyakov loop in chiral quark models at finite temperature, Phys. Rev. D 74, 065005 (2006), arXiv:hep-ph/0412308 .
  * Lo  _et al._ [2021] P. M. Lo, K. Redlich, and C. Sasaki, Fluctuations of the order parameter in an $SU(N_{c})$ effective model, Phys. Rev. D 103, 074026 (2021), arXiv:2101.12663 [hep-ph] .
  * Mocsy  _et al._ [2004] A. Mocsy, F. Sannino, and K. Tuominen, Confinement versus chiral symmetry, Phys. Rev. Lett. 92, 182302 (2004), arXiv:hep-ph/0308135 .
  * Ratti  _et al._ [2006] C. Ratti, M. A. Thaler, and W. Weise, Phases of QCD: Lattice thermodynamics and a field theoretical model, Phys. Rev. D 73, 014019 (2006), arXiv:hep-ph/0506234 .
  * Hansen  _et al._ [2007] H. Hansen, W. M. Alberico, A. Beraudo, A. Molinari, M. Nardi, and C. Ratti, Mesonic correlation functions at finite temperature and density in the Nambu-Jona-Lasinio model with a Polyakov loop, Phys. Rev. D 75, 065004 (2007), arXiv:hep-ph/0609116 .
  * Moreira  _et al._ [2018] J. Moreira, J. Morais, B. Hiller, A. A. Osipov, and A. H. Blin, Thermodynamical properties of strongly interacting matter in a model with explicit chiral symmetry breaking interactions, Phys. Rev. D 98, 074010 (2018), arXiv:1806.00327 [hep-ph] .
  * Câmara Pereira  _et al._ [2020] R. Câmara Pereira, J. a. Moreira, and P. Costa, The strange critical endpoint and isentropic trajectories in an extended PNJL model with eight Quark interactions, Eur. Phys. J. A 56, 214 (2020), arXiv:2006.02385 [hep-ph] .
  * Pereira [2021] R. C. Pereira,  _Quantum Chromodynamics Phase Diagram Under Extreme Conditions_ , Ph.D. thesis, Coimbra U. (2021).
  * Ferreira  _et al._ [2021] M. Ferreira, R. Câmara Pereira, and C. Providência, Hybrid stars with large strange quark cores constrained by GW170817, Phys. Rev. D 103, 123020 (2021), arXiv:2105.06239 [nucl-th] .
  * Haas  _et al._ [2013] L. M. Haas, R. Stiele, J. Braun, J. M. Pawlowski, and J. Schaffner-Bielich, Improved Polyakov-loop potential for effective models from functional calculations, Phys. Rev. D 87, 076004 (2013), arXiv:1302.1993 [hep-ph] .
  * Lo  _et al._ [2018] P. M. Lo, M. Szymański, K. Redlich, and C. Sasaki, Polyakov loop fluctuations in the presence of external fields, Phys. Rev. D 97, 114006 (2018), arXiv:1801.08040 [hep-ph] .
  * Rößner [2006] S. Rößner,  _Field theoretical modelling of the QCD phase diagram_ , Ph.D. thesis, Technische Universität München (2006).
  * Hell  _et al._ [2010] T. Hell, S. Rossner, M. Cristoforetti, and W. Weise, Thermodynamics of a three-flavor nonlocal Polyakov-Nambu-Jona-Lasinio model, Phys. Rev. D 81, 074034 (2010), arXiv:0911.3510 [hep-ph] .
  * Hell [2010] T. Hell,  _Modeling the Thermodynamics of QCD_ , Ph.D. thesis, Technische Universität München (2010).
  * Sasaki and Redlich [2012] C. Sasaki and K. Redlich, An Effective gluon potential and hybrid approach to Yang-Mills thermodynamics, Phys. Rev. D 86, 014007 (2012), arXiv:1204.4330 [hep-ph] .
  * Hell  _et al._ [2009] T. Hell, S. Roessner, M. Cristoforetti, and W. Weise, Dynamics and thermodynamics of a non-local PNJL model with running coupling, Phys. Rev. D 79, 014022 (2009), arXiv:0810.1099 [hep-ph] .
  * Drouffe and Zuber [1983] J.-M. Drouffe and J.-B. Zuber, Strong coupling and mean field methods in lattice gauge theories, Physics Reports 102, 1 (1983).