# Metastability and topology in the magnetic topological insulator MnBi2Te4

Jeonghwan Ahn,1,∗, Seoung-Hun Kang1,∗, Mina Yoon1, Panchapakesan Ganesh,2,‡
Jaron T. Krogel1,† Materials Science and Technology Division, Oak Ridge
National Laboratory, Oak Ridge, Tennessee 37831, USA  Center for Nanophase
Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831,
USA

###### Abstract

We study the effect of stacking faults on the topological properties of the
magnetic topological insulator MnBi2Te4 (MBT) using density functional theory
calculations and the Hubbard $U$ being tuned with many-body diffusion Monte
Carlo techniques. We show that a modest deviation from the equilibrium
interlayer distance leads to a topological phase transition from a non-trivial
to a trivial topology, suggesting that tuning the interlayer coupling by
adjusting the interlayer distance alone can lead to different topological
phases. Interestingly, due to the locally increased interlayer distance of the
top layer, a metastable stacking fault in MBT leads to a nearly gapless state
at the topmost layer due to charge redistribution as the topmost layer
recedes. We further find evidence of spin-momentum locking in the surface
state along with a weak preservation of the band inversion in the near gapless
state, which is indicative of the non-trivial topological surface states for
the metastable stacking fault. Our findings provide a possible explanation for
reconciling the long-standing puzzle of gapped and gapless states on MBT
surfaces.

∗ These authors contributed equally to this work.

†<EMAIL_ADDRESS>‡<EMAIL_ADDRESS>
## Introduction

MnBi2Te4 (MBT) has drawn a great deal of interest as the first realization of
an intrinsic and stoichiometric topological insulator with $A$-type
antiferromagnetic (AFM) order Otrokov et al. (2019a). The structure of MBT is
characterized by a layered structure consisting of septuple layers coupled
together with van der Waals (vdW) interactions as illustrated in Fig. 1(a).
Because of its crystal structure, there is a topological invariant under the
product operation $S$ of time reversal ($\Theta$) and primitive lattice
translation ($T_{1/2}$), namely $S=\Theta T_{1/2}$, giving rise to $Z_{2}$
topological classification Mong et al. (2010). Consequently, gapped surface
states are predicted to exist at the $S$-breaking surface of (0001), which
makes MBT a highly desired platform as a single material to study an interplay
between magnetism and topology. Broken time-reversal symmetry in MBT leads to
novel topological phases such as quantum anomalous Hall effect Deng et al.
(2020); Zhao et al. (2021) and axion insulator Liu et al. (2020); Lüpke et al.
(2022).

However, the existence of the gapped surface state has remained elusive in the
sense that there has been no clear consensus among a series of experimental
studies; the measured surface gaps differ depending on the sample Shikin et
al. (2021); Garnica et al. (2022) and even gapless surface states have been
reported Hao et al. (2019); Chen et al. (2019a); Li et al. (2019a); Swatek et
al. (2020); Nevola et al. (2020) since the first reports that support the
gapped surface states of MBT Otrokov et al. (2019a). Furthermore, a sizable
spread among the measured surface gaps brought a fundamental question
regarding whether the surface gap originated from the AFM order or not Chen et
al. (2019a). There have been many attempts to elucidate observed gapless and
sample-dependent topological surface states, among which the different surface
magnetic orderings from bulk ones like G-type or in-plane A-type AFM are
predicted to result in linear dispersion Hao et al. (2019); Chen et al.
(2019a). Defects have been considered one of the probable sources to induce
the surface magnetism distinct from the bulk one because of frequent
observation of various defects in many experimental studies Zeugner et al.
(2019); Huang et al. (2020); Lai et al. (2021); Garnica et al. (2022).
Especially, the Mn-involved defects are known to result in ferrimagnetic
structure in the septuple layers Lai et al. (2021), which are understood to
shift the real-space localization of the topological surface states compared
to the ideal case and lead to a substantial reduction of the surface gap
Garnica et al. (2022). Besides, the topological properties were sensitive to
the hydrostatic pressure (strain) Chen et al. (2019b); Guo et al. (2021); Xu
et al. (2022) because structural and magnetic properties are found to be
highly coupled with pressure, reflecting their combined role in determining
the topology of the system. Nevertheless, none of the above fully resolves the
ongoing debate on surface gaps.

For Te-based layered materials, strong interlayer coupling due to interlayer
hybridization is known to host correlated physics such as the metal-insulator
transition between bilayer and monolayer IrTe2 with the formation of charge
density wave phase for monolayer Hwang et al. (2022), and Lifshitz transition
for PtTe2 Lin et al. (2020). One can also expect substantial interlayer
coupling in the MBT system because its valence band maximum (VBM) also mostly
consists of pz orbitals. Hence interlayer hybridization can be maximized
between pz orbitals aligned with the interlayer axis. More interestingly,
since not only Te-pz, but also Bi-pz orbitals are involved in the VBM of MBT,
the spin-orbit coupling (SOC) effects would play some role in the interlayer
hybridization between adjacent Te planes, which has not been the focus of the
previous literature. The SOC-involved interlayer coupling would directly
affect the band inversion and, thus, the topology of the system. Furthermore,
the fact that the topology of MBT becomes nontrivial as the layers are
stacked, but the single-layer MBT is topologically trivial Otrokov et al.
(2019b) underscores the interlayer coupling as a key factor to drive the non-
trivial topology of the system.

We note that the interlayer coupling can be tuned by the nature of the local
stacking, such as different stacking modes and twist angles, because of
different resulting interlayer separations and the nearest neighbors across
the vdW gap, which results in different physical properties. Indeed, according
to a recent experiment Yan (2022), stacking faults induced by repeated bending
of MBT bulk samples are observed to suppress the Néel temperature compared to
the ideal case. The suppression effect receded after about one hour, is
indicative of the metastability of the stacking faults induced by the bending.
Furthermore, it is remarkable to note that the impact of the stacking fault
(or metastable stacking configuration) on the electronic structure, especially
the topology, in MBT has not been explored yet though its existence is
strongly suggestive in the mother material Bi2Te3 Medlin et al. (2019) or MBT
variant such as MnBi4Te7 or MnBi6Te10 Yan et al. (2020).

The above facts motivated us to explore in greater depth the possible
existence of a metastable stacking order in MBT and its influence on surface
gap formation and corresponding topology. We firstly show that SOC and vdW
interactions are intertwined due to interlayer hybridization involved with Bi
pz and Te -pz, supported by different predictions of interlayer binding
energies and separations depending on the inclusion of SOC. Next, we
demonstrate that modest expansion of the interlayer separations results in a
topological phase transition, indicating the significance of the interlayer
coupling to affect the topological phase, even with modest perturbations to
the lattice. Furthermore, we show that a metastable local stacking mode in MBT
characterized by a local increase in interlayer separation that lies right at
the verge of creating inverted bands. As demonstrated in a 6-layer MBT system,
we find that a metastable stacking at the MBT surface (which can be formed by
sliding the topmost layer) leads to nearly gapless states localized at the
topmost layer, consistent with the experimental observation of the gapless
surface state for bulk MBT Lai et al. (2021); Garnica et al. (2022). The spin
is found to be locked to the momentum at the surface state and the gap closure
at the surface corresponds to band inversion between the surface CBM and VBM,
which is suggestive of its non-trivial topological surface states.

## Results

### Interplay between SOC and vdW interactions on interlayer hybridization

Figure 1: (a) Crystal structure of MBT. The dotted lines represent its unit
cell. (b) Interlayer binding energies for the bulk MBT as a function of the
vertical distance between adjacent Te layers with/without SOC. The black solid
(open) symbol represents interlayer binding energy with (without) SOC. (c) The
Bandgap at $\Gamma$ of bulk MBT depends on including SOC as a function of the
vertical distance between adjacent Te layers. The black open symbol represents
the data without SOC, while the red solid (open) symbol denotes SOC data for
which the corresponding topological $Z_{2}$ index is 1 (0). Band structures of
the bulk MBT with the inclusion of SOC at (d) the equilibrium interlayer
separation of 2.74 Å and (e) 3.4 Å.

We first discuss the role of the SOC interaction in determining the structural
and interlayer binding properties of the MBT layers. To this end, we present
Fig. 1(b) to display interlayer binding energies of bulk MBT as a function of
the vertical distance between adjacent Te layers, interlayer separation
($d_{\text{Te-Te}}$), with and without the inclusion of SOC at PBE+D3 level.
The binding curves are found to be distinct depending on the inclusion of SOC.
The equilibrium binding distance and interlayer separation for the SOC
calculations are estimated to be 2.74 Å and 21.3 meV/Å2, respectively. The
predicted interlayer separation from the SOC calculations agrees with
experimental values based on an X-ray measurement of 2.72(1) Å Zeugner et al.
(2019). On the other hand, for the non-SOC calculations, the corresponding
quantities are 2.92 Å and 18.94 meV/Å2, which are smaller (longer) than the
non-SOC binding energy (interlayer separation) about 11 $\%$ (7 $\%$). This
discrepancy depending on the inclusion of SOC, indicates that SOC and vdW
interactions are intertwined, which has not been the focus in previous
literature. The intertwining of SOC and vdW interactions is found to be
regardless of the inclusion of Hubbard U, the value of U, and different
pseudopotentials.

In comparison with the interlayer binding curves of materials in a similar
class as MBT, as shown in Supplementary Fig.1, the interlayer binding curves
for the Bi-involved compounds display the difference between SOC and non-SOC
calculations while the binding curves are identical for the Sb-involved
compounds. This suggests that the Bi element possessing larger SOC strength
than Sb is responsible for the intertwining SOC and vdW interactions. Mn does
not play a significant role in the intertwining effect since binding curves
for the Mn-involved compound show analogous features to those for non-magnetic
counterparts. The difference in chalcogen atoms is sufficient to explain the
different curvature of the interlayer binding curves.

To further understand the interplay between SOC and vdW interactions, we
investigated band structures and the band gap evolution as a function of
$d_{\text{Te-Te}}$. Figure 1(d) displays band structures of the bulk MBT at
the equilibrium interlayer separation when including SOC. Both orbital
contributions of mainly Bi-pz and Te-pz to VBM signify the band inversion.
Considering that this interlayer hybridization is aligned with the axis along
which MBT septuple layers are stacked through the vdW interactions, the
increasing SOC-involved interlayer hybridization as the septuple layers become
close to each other is understood to result in the enhanced interlayer
bindings in MBT when incorporating SOC. A trend of band gap as a function of
interlayer separation shown in Fig. 1(c) is considered another clue indicating
the coupling between SOC and vdW interactions in the interlayer hybridization.
Figure 1(c) displays band gaps at $\Gamma$ point for different interlayer
separations depending on the inclusion of SOC. When SOC is not included, the
band gap is observed to be monotonically reduced with decreasing interlayer
distances, while the band gap decreases as interlayer separation decreases but
increases at the region of $d_{\text{Te-Te}}<3.2$ in the presence of SOC. We
attribute the SOC band gap trend distinct from the non-SOC case to the
topological phase transition for the interlayer separations, which is verified
with the $Z_{2}$ topological invariant changed from 1 to 0 between 3.1 $\sim$
3.2 Å as shown in Fig. 1(c). This is in line with the bulk gap closure between
3.1 $\sim$ 3.2 Å and the band inversion disappears with increasing interlayer
distance as shown in Fig. 1(e). Interestingly, we notice that the band gap
trend with SOC is identical to the band gap trend as a function of the SOC
strength (See Fig.2D of Ref. Li et al. (2019b)) where phase transition occurs
from the topologically non-trivial state to normal insulator. This suggests
that interlayer separation controls the SOC strength through interlayer
hybridization, and the modest deviation from its equilibrium value could
result in a significant change in topological properties. From the above
facts, we conclude that there is huge sensitivity in topological properties to
the interlayer separation, and it is one of the most crucial factors
responsible for the topological properties of MBT.

### Metastable stacking at the MBT bulk surface

Figure 2: (a) Potential energy surfaces for trilayer MBT with respect to
sliding its topmost septuple layer with the inclusion of SOC. The letters A,
B, C, and M denote the stacking position of the topmost septuple layer on top
of the AB-stacked bilayer. (b) The DFT total energies as a function of
interlayer separation between the adjacent Te plane of the topmost and second-
topmost septuple layers, $d_{\text{Te-Te}}$, relative to the energy minima for
the ABC- and ABCABC-stacked layers, respectively. Third-order polynomial fits
are presented for a guide to the eye.

Considering that the topological properties of MBT might be tuned by stacking
modes that have different interlayer separations from the ideal case based on
the above results, we investigated the possible existence of metastable
stacking mode for bulk MBT. To this end, we first computed the potential
energy surface (PES) for trilayer MBT with respect to the sliding of its
topmost layer on top of the underlying AB-stacked bilayer as displayed in Fig.
2(a). For the consistency of the above binding curves for bulk MBT, the PES
was computed with PBE+$U$(3.5 eV)+D3+SOC to reflect the intertwining of vdW
and SOC interactions utilizing total energy curves as a function of interlayer
separation between the topmost and second topmost layers at the given sliding
with the thickness of MBT being fixed with the bulk value as shown in Fig.
2(a). The letters A, B, and C in the contour map represent the relative
stacking of the topmost layer above the AB-stacked one, and hence C appears to
be the lowest-energy point because the ground-state stacking for the bulk
structure is ABC. As shown in Fig. 2(b), we could also identify the ABA
stacking mode as a local minimum in the energy map with the energy difference
from the ground-state structure being $\sim$ 2 meV/ Å2, which corresponds to
about 11 $\%$ of the interlayer binding energy of the bulk MBT. We find from
third-order polynomial fits that the equilibrium interlayer separation at ABA
stacking has increased by about 11 $\%$ with an energy minimum residing at
$d_{\text{Te-Te}}$ = 3.12 Å, relative to ABC stacking whose minimum is located
at $d_{\text{Te-Te}}$ = 2.74 Å while the transition state (ABM stacking)
exhibits the energy minimum at $d_{\text{Te-Te}}$ = 3.18 Å. Furthermore, one
can clearly see in the PES energetically favored transition path from the
global minimum to the local one (See path C-M-A Fig. 2(a)) along which the
corresponding barrier is estimated to be $\sim$ 3 meV/ Å2. This barrier and
the small energy difference between global and local minima indicate that the
metastable stacking mode might easily occur.

Noting that 6-layer MBT is regarded as a sufficiently large slab model to
approximate electronic structures of the surface of bulk MBT in the previous
literature Shikin et al. (2021), we examined with the 6-layer MBT whether
metastable stacking of the topmost layer could occur at the surface of bulk
MBT. Provided that the transition path for the 6-layer MBT, whose topmost
layer is on top of the ABCAB-stacked one, is identical to that of the trilayer
MBT, we compared their total energies relative to the lowest-energy stacking
ones as a function of interlayer separation between the topmost and second
topmost layers, $d_{\text{Te-Te}}$ when the topmost one is located at A, B or
C as shown in Fig. 2(b). The total energy curves of trilayer MBT are confirmed
to be nearly identical to that of 6-layer MBT for all slidings, and
consequently, energy differences between the stacking modes and corresponding
interlayer separations become nearly identical to each other, with the
equilibrium interlayer separation of the topmost layer residing at
$d_{\text{Te-Te}}$ = 3.12 Å. We also point out that slab and bulk calculations
yield the same interlayer separations of 2.74 Å for the lowest energy ABC and
ABCABC stacking mode, suggesting that the interlayer separation is not
sensitive to the number of MBT septuple layers because of its substantial
thickness ($\sim 10.74$ Å). This justifies our approach with a 6-layer MBT for
energetic investigation of the possible metastable stacking mode at the bulk
surface. From this, we predict that these equilibrium quantities, along with
the transition path to the metastable stacking mode, would be preserved with
increasing thickness of MBT so that metastable stacking mode could also easily
occur at the surface of bulk MBT with a transition barrier of $\sim$ 3 meV/
Å2.

Figure 3: Band structures for (a) the ground-state (ABCABC) and (d) metastable
stacking (ABCABA) modes at their equilibrium. The solid violet symbols
represent the weight of the projections of the Bloch states on the atoms
belonging to the topmost and bottommost septuple layers. Spatial distributions
of the CBM and VBM states at the gamma point are represented in (b) and (e)
for the ground-state and the metastable stacking modes, respectively. The
isosurfaces levels are $\pm 0.001$ $e$/bohr3 and only two topmost layers are
presented for clarity. Spin maps plotted in the reciprocal space for the CBM
state (c) for the ground-state stacking and (f) metastable stacking modes. The
size of black arrows and different colors represent the in-plane and out-of-
plane spin components.

### Electronic structure of metastable stacking mode

We now present the band structures for the ground-state and metastable
stacking modes identified in Fig. 2. Fig. 3(a) presents the orbital
projections of the surface atoms for the ground-state stacking at the
equilibrium separation, defined here as those residing in the top- and bottom-
most septuple layers, out of the total band structures denoted by violet
symbols. It is found that occupied levels with strong surface contributions
generally lie below sub-surface dominant states that are largely bulk-like,
indicating that CBM and VBM originate from surface and bulk-like parts,
respectively, as displayed in Fig. 3(b) showing the spatial distributions of
CBM and VBM states at the gamma point. The surface gap of 118 meV is close to
the experimentally measured gaps ranging from 20 $\sim$ 100 meV Otrokov et al.
(2019a); Lee et al. (2019); Vidal et al. (2019); Zeugner et al. (2019); Shikin
et al. (2021); Garnica et al. (2022) while the bulk gap of 205 meV is found to
be quantitatively in agreement with the experimental values ranging from 180
$\sim$ 220 meV Chen et al. (2019c). Furthermore, we present in Fig. 3(c) the
spin projection onto the lateral plane of the reciprocal space for the
surface-like CBM state, showing the spin is largely locked to the momentum, as
is commonly observed in topologically non-trivial surface states. This
qualitative observation is also consistent with the non-trivial topology of
the bulk MBT at the equilibrium distance confirmed in Fig. 1(c).

On the other hand, the band structure of the metastable stacking mode is seen
to be noticeably different compared with the ground-state one, as shown in
Fig. 3(d); the surface-atom contributions appear in VBM as a result of the
expanded separation distance of the topmost layer while CBM still remains as a
surface band (See also Fig. 3(e)). Consequently, the surface gap is estimated
to be 5 meV at the metastable stacking mode, which may explain the occasional
experimental findings of the gapless surface state Hao et al. (2019); Chen et
al. (2019a). We note that CBM for metastable stacking can still be seen as the
surface-like state and the corresponding orbital character of is qualitatively
analogous to that of the ground-state stacking whose surface state is in
topologically non-trivial (See Fig. 3(b) and (e)) in the sense that the CBM
state is mainly of Te-pz and Bi-pz. Furthermore, we also find spin-momentum
locking for the metastable stacking mode, as shown in Fig. 3(f) through the
corresponding signal is observed to be weakened compared to the case of
ground-state stacking. This suggests that the near gapless metastable stacking
might still host a topologically non-trivial surface state even in the
presence of weaker interlayer coupling.

Figure 4: Band structures of the ABCABC-stacked MBT layer for different
$d_{\text{Te-Te}}$, that is, (a) $d_{\text{Te-Te}}$ = 2.74 Å, (b)
$d_{\text{Te-Te}}$ = 2.94 Å, (c) $d_{\text{Te-Te}}$ = 3.0 Å and (d)
$d_{\text{Te-Te}}$ = 3.1 Å. The solid violet symbols represent the weight of
the projections of the Bloch states on the atoms belonging to the topmost and
bottommost septuple layers. The insets represent Spatial distributions of the
CBM and VBM states at the gamma point at each separation. The isosurfaces
levels are $\pm 0.001$ $e$/bohr3 and only two topmost layers are presented for
clarity.

Since the significant change in the band structure at the metastable stacking
mode is understood to be due to the change of interlayer coupling by the
expansion of the equilibrium separation at the topmost layer, we explore the
effect of interlayer separation alone for the topmost layer on the band
structure of the 6-layer MBT slab. To this end, we analyze in Fig. 4 an
evolution of the band structures of the ground-state stacking mode by varying
only interlayer separation of the topmost layer, $d_{\text{Te-Te}}$. As the
topmost layer recedes, energy levels near the Fermi level significantly
change. Especially, both VBM and CBM move toward the center of the gap, with
the VBM character changing from bulk-like to surface-dominant while the CBM
maintains strong surface contributions, giving rise to a significant band gap
reduction of the surface state. Consequently, when the top-most layer reaches
a separation of $d_{\text{Te-Te}}$ = 3.1 Å (See Fig. 4(e)), both CBM and VBM
mainly have the surface contributions and the surface gap is estimated to be
11 meV as similar to the case of the metastable stacking configuration in Fig.
3(d). This indicates that the surface states (gap) are very sensitive to the
interlayer separation between the topmost and second-topmost layers, which
could be considered a crucial factor in explaining gapless or sample-dependent
surface gap measured in the previous experiments Shikin et al. (2021); Garnica
et al. (2022); Hao et al. (2019); Chen et al. (2019a).

As is clearly seen in the spatial distributions of CBM and VBM at the gamma
point presented in the inset of Fig. 4, the CBM state remains localized at the
interface of the topmost layer, while the localization of the VBM state
changes from the bulk to the topmost region as the topmost layer recedes. Thus
the gapless state, if any, is expected to exist primarily at the topmost layer
exposed to the spectroscopic measurement’s light source. Furthermore, we
notice the band inversion occurring between $d_{\text{Te-Te}}$ = 3.0 and 3.1 Å
since the spatial distributions of CBM and VBM are observed to be switched.
Analysis of the orbital projections reveals that VBM and CBM mainly have Te-pz
and Bi-pz character, respectively, at $d_{\text{Te-Te}}$ = 3.0 and 3.1 Å and
these are seen to be inverted at $d_{\text{Te-Te}}$ = 3.1 Å. The separation
range of the disappearance of band inversion, which coincides with the one for
the band-gap closure, might suggest a topological phase transition occurring
between $d_{\text{Te-Te}}$ = 3.0 and 3.1 Å. This requires band-gap closure at
the intermediate critical point, which is consistent with the gapless surface
state Lai et al. (2021); Garnica et al. (2022). We further see that the
splitting in the VBM states is about 65 meV near $d_{\text{Te-Te}}$ = 3.1 Å in
the ABCABC stacking, reproducing well the observed 70 meV $k_{z}$ dispersion
in the VBM states for the MBT system with a gapless surface state Hao et al.
(2019). In addition, it is worthwhile to note that the total band structures
of ABCABA and ABCABC at $d_{\text{Te-Te}}$ = 3.1 Å are found to be nearly
identical to each other, but the band is inverted at the topmost layer for the
ABCABA stacking while there is no band inversion for the ABCABC one (See Fig.
3(e) and Fig. 4(d)). This suggests that changes in interlayer hybridization in
the vicinity of $d_{\text{Te-Te}}$ = 3.1 Å for the two stackings are
sufficient to flip the band inversion when it is close to forming; however,
the gap closing is clearly dominated by changes in interlayer coupling driven
by the separation distance, underscoring the significance of the interlayer
coupling in understanding the surface gap closure and its corresponding
topological phase.

## Discussion

Combined, our results suggest that local realizations of this metastable
stacking mode in MBT result in near-gap closings, which leads us to conclude
that the proposed metastable stacking mode could explain the controversial
spread in the experimentally measured surface gap. It is worthwhile to note
that tuning interlayer separation was recently achieved through an
electrochemical protocol in MoS2 by introducing active sites of Ce ions and
intercalating Pt atoms in the interlayer region Ding et al. (2021). This makes
us expect the interlayer separation of MBT to be also tuned in a similar
fashion, which opens a possibility to manipulate its topological phase via
direct control of interlayer separation.

We have proposed the metastable stacking mode that could consistently explain
the spread in the experimentally measured surface gap with the non-trivial
topology. Especially in light of our results, the surface-gap closure may be
understood as a part of a topological phase transition induced by the
expansion of the interlayer separation at a metastable stacking configuration.
The interplay between the interlayer hybridization of the pz orbitals of Bi
and Te and the effective spin-orbit coupling strength is important to this
effect. This reveals the crucial role of interlayer separation as the tuning
factor of interlayer coupling, which directly impacts the energetics, band
structure, and the topology of MBT. In this sense, we point out that increased
interlayer separation could be induced by surface defects Zeugner et al.
(2019); Huang et al. (2020); Lai et al. (2021); Garnica et al. (2022) or
intercalation Mazza et al. (2022), which were proposed to elucidate the
gapless surface state in the previous literature. Thus, our study provides
deeper insight into understanding the interplay between structure and topology
while shedding light on the possibility of manipulating the topological phase
via engineering the interlayer separations.

## Methods

The DFT calculations were implemented with the QUANTUM ESPRESSO package
Giannozzi et al. (2009). To solve the Kohn-Sham equations, we used optimized
norm-conserving Vanderbilt (ONCV) pseudopotentials Hamann (2013) for the
atomic species in MBT. The plane-wave cutoff was set to be 120 Ry, and the
$15\times 15\times 3$ ($15\times 15\times 1$) Monkhorst–Pack Monkhorst and
Pack (1976) $k$ meshes were used for the bulk (few-layer) MBT. The structural
relaxations were done with the accuracy of 10-4 Ry/Bohr. The van der Waals
interactions were included with DFT+D3 method Grimme et al. (2010), and the
Hubbard U Dudarev et al. (1998); Liechtenstein et al. (1995) was applied to
properly describe the Mn-driven magnetic states. The spin-orbit interactions
were also included to discuss the effects of spin-orbit coupling (SOC) on the
electronic and structural properties and the energetics at the DFT level. In
structural models of few-layer MBT systems, we used a vacuum distance of 60 Å
along the vertical direction up to the trilayer while 120 Å is used for
6-layer MBT to minimize spurious interactions between periodic images. From
the DFT results, we obtain the maximally localized Wannier functions for
4$s$\- and 3$d$-orbitals of Mn atom, 6$p$ of Bi atom and 5$p$-orbitals on Te
atoms by using the WANNIER90 code Mostofi et al. (2014). These were used to
analyze the topological properties. The topological index are calculated by
the WANNIERTOOLS Wu et al. (2018), which is based on the Wannier charge center
calculation method Soluyanov and Vanderbilt (2011).

We utilized QMCPACK code for the fixed-node DMC calculations to tune Hubbard U
for the monolayer MBT as in our previous study for the bulk MBT Bennett et al.
(2022). The pseudopotentials which were generated for Mn Annaberdiyev et al.
(2018) or tested for Bi Metz et al. (2000) and Te Stoll et al. (2002) within
the ccECP scheme Bennett et al. (2017); Annaberdiyev et al. (2018), yielded
bulk Mn magnetization in benchmark quality. Slater-Jastrow type trial
wavefunctions were used with the Jastrow parameters up to three-body ones to
incorporate electron-ion, electron-electron, and electron-electron-ion
correlations, which were optimized with variational Monte Carlo calculations
based on the linear method of Umrigar et al. Umrigar et al. (2007). In
subsequent DMC calculations, we used a time step of $\tau=0.005$ Ha-1 and
size-consistent T-moves for variational evaluation of non-local
pseudopotentials when imaginary time projection proceeds Kim et al. (2018);
Kent et al. (2020).

To optimize the nodal surface of the many-body wavefunctions for the monolayer
MBT, we computed its DMC total energies for the 7-atom ferromagnetic cell as a
function of U, for which the PBE+U orbitals are used in constructing the
wavefunctions. Fitting DMC energies to the third-order polynomial fit results
in the optimal value of U= 3.7(3) eV (See Supplementary Fig.2). The DMC-
optimized U for the monolayer is observed to be identical to the corresponding
value for the bulk MBT within the statistical error bar Bennett et al. (2022),
which provides the onsite Mn magnetization that is in good agreement with
neutron diffraction measurement Ding et al. (2020). The nearly identical U
values for bulk and monolayer MBT are understood by the similar chemical
environment around Mn atoms, even though layer thickness varies. Based on this
fact, the following DFT calculations for few-layer MBT were performed with the
DMC-informed Hubbard U for the bulk value of 3.5 eV for the consistency of our
calculations dealing with both few-layer and bulk MBT.

## Acknowledgments

Work performed by J.A., P.G., and J.T.K. (DFT calculations, original idea) was
supported by the U.S. Department of Energy, Office of Science, Basic Energy
Sciences, Materials Sciences and Engineering Division, as part of the
Computational Materials Sciences Program and Center for Predictive Simulation
of Functional Materials. Work performed by S.-H. K. and M. Y. (topological
index and spin texture calculations) was supported by the U.S. Department of
Energy (DOE), Office of Science, National Quantum Information Science Research
Centers, Quantum Science Center. An award of computer time was provided by the
Innovative and Novel Computational Impact on Theory and Experiment (INCITE)
program. This research used resources of the Oak Ridge Leadership Computing
Facility, which is a DOE Office of Science User Facility supported under
Contract DE-AC05-00OR22725. This research also utilized resources of the
National Energy Research Scientific Computing Center, a DOE Office of Science
User Facility supported by the Office of Science of the U.S. Department of
Energy under Contract No. DE-AC02-05CH11231 using NERSC award BES-
ERCAP0024568.

## Data Availability

The data that support the findings of this study are available within this
article, its supplementary material, and in the materials data
facilityBlaiszik et al. (2016, 2019) at [link to be provided upon acceptance
of this manuscript].

## Competing interests

The authors declare no competing interests.

## References

## References

  * Otrokov et al. (2019a) M. M. Otrokov, I. I. Klimovskikh, H. Bentmann, D. Estyunin, A. Zeugner, Z. S. Aliev, S. Gaß, A. Wolter, A. Koroleva, A. M. Shikin, et al., Nature 576, 416 (2019a).
  * Mong et al. (2010) R. S. Mong, A. M. Essin, and J. E. Moore, Phys. Rev. B 81, 245209 (2010).
  * Deng et al. (2020) Y. Deng, Y. Yu, M. Z. Shi, Z. Guo, Z. Xu, J. Wang, X. H. Chen, and Y. Zhang, Science 367, 895 (2020).
  * Zhao et al. (2021) Y.-F. Zhao, L.-J. Zhou, F. Wang, G. Wang, T. Song, D. Ovchinnikov, H. Yi, R. Mei, K. Wang, M. H. Chan, et al., Nano Lett. 21, 7691 (2021).
  * Liu et al. (2020) C. Liu, Y. Wang, H. Li, Y. Wu, Y. Li, J. Li, K. He, Y. Xu, J. Zhang, and Y. Wang, Nat. Mater. 19, 522 (2020).
  * Lüpke et al. (2022) F. Lüpke, A. D. Pham, Y.-F. Zhao, L.-J. Zhou, W. Lu, E. Briggs, J. Bernholc, M. Kolmer, J. Teeter, W. Ko, et al., Phys. Rev. B 105, 035423 (2022).
  * Shikin et al. (2021) A. M. Shikin, D. Estyunin, N. L. Zaitsev, D. Glazkova, I. I. Klimovskikh, S. Filnov, A. G. Rybkin, E. Schwier, S. Kumar, A. Kimura, et al., Phys. Rev. B 104, 115168 (2021).
  * Garnica et al. (2022) M. Garnica, M. M. Otrokov, P. C. Aguilar, I. Klimovskikh, D. Estyunin, Z. S. Aliev, I. R. Amiraslanov, N. A. Abdullayev, V. N. Zverev, M. B. Babanly, et al., npj Quantum Mater. 7, 1 (2022).
  * Hao et al. (2019) Y.-J. Hao, P. Liu, Y. Feng, X.-M. Ma, E. F. Schwier, M. Arita, S. Kumar, C. Hu, R. Lu, M. Zeng, et al., Phys. Rev. X 9, 041038 (2019).
  * Chen et al. (2019a) Y. J. Chen, L. X. Xu, J. H. Li, Y. W. Li, H. Y. Wang, C. F. Zhang, H. Li, Y. Wu, A. J. Liang, C. Chen, et al., Phys. Rev. X 9, 041040 (2019a).
  * Li et al. (2019a) H. Li, S.-Y. Gao, S.-F. Duan, Y.-F. Xu, K.-J. Zhu, S.-J. Tian, J.-C. Gao, W.-H. Fan, Z.-C. Rao, J.-R. Huang, et al., Phys. Rev. X 9, 041039 (2019a).
  * Swatek et al. (2020) P. Swatek, Y. Wu, L.-L. Wang, K. Lee, B. Schrunk, J. Yan, and A. Kaminski, Phys. Rev. B 101, 161109 (2020).
  * Nevola et al. (2020) D. Nevola, H. X. Li, J.-Q. Yan, R. Moore, H.-N. Lee, H. Miao, and P. D. Johnson, Phys. Rev. Lett. 125, 117205 (2020).
  * Zeugner et al. (2019) A. Zeugner, F. Nietschke, A. U. Wolter, S. Gaß, R. C. Vidal, T. R. Peixoto, D. Pohl, C. Damm, A. Lubk, R. Hentrich, et al., Chem. Mater. 31, 2795 (2019).
  * Huang et al. (2020) Z. Huang, M.-H. Du, J. Yan, and W. Wu, Phys. Rev. Mater. 4, 121202(R) (2020).
  * Lai et al. (2021) Y. Lai, L. Ke, J. Yan, R. D. McDonald, and R. J. McQueeney, Phys. Rev. B 103, 184429 (2021).
  * Chen et al. (2019b) K. Chen, B. Wang, J.-Q. Yan, D. Parker, J.-S. Zhou, Y. Uwatoko, and J.-G. Cheng, Phys. Rev. Mater. 3, 094201 (2019b).
  * Guo et al. (2021) W.-T. Guo, L. Huang, Y. Yang, Z. Huang, and J.-M. Zhang, New J. Phys. 23, 083030 (2021).
  * Xu et al. (2022) Z. Xu, M. Ye, J. Li, W. Duan, and Y. Xu, Phys. Rev. B 105, 085129 (2022).
  * Hwang et al. (2022) J. Hwang, K. Kim, C. Zhang, T. Zhu, C. Herbig, S. Kim, B. Kim, Y. Zhong, M. Salah, M. M. El-Desoky, et al., Nat. Commun. 13, 1 (2022).
  * Lin et al. (2020) M.-K. Lin, R. A. B. Villaos, J. A. Hlevyack, P. Chen, R.-Y. Liu, C.-H. Hsu, J. Avila, S.-K. Mo, F.-C. Chuang, and T.-C. Chiang, Phys. Rev. Lett. 124, 036402 (2020).
  * Otrokov et al. (2019b) M. Otrokov, I. P. Rusinov, M. Blanco-Rey, M. Hoffmann, A. Y. Vyazovskaya, S. Eremeev, A. Ernst, P. M. Echenique, A. Arnau, and E. V. Chulkov, Phys. Rev. Lett. 122, 107202 (2019b).
  * Yan (2022) J.-Q. Yan, ECS J. Solid State Sci. Technol. 11, 063007 (2022).
  * Medlin et al. (2019) D. L. Medlin, N. Yang, C. D. Spataru, L. Hale, and Y. Mishin, Nat. Commun. 10, 1 (2019).
  * Yan et al. (2020) J.-Q. Yan, Y. Liu, D. Parker, Y. Wu, A. Aczel, M. Matsuda, M. McGuire, and B. Sales, Phys. Rev. Mater. 4, 054202 (2020).
  * Li et al. (2019b) J. Li, Y. Li, S. Du, Z. Wang, B.-L. Gu, S.-C. Zhang, K. He, W. Duan, and Y. Xu, Sci. Adv. 5, eaaw5685 (2019b).
  * Lee et al. (2019) S. H. Lee, Y. Zhu, Y. Wang, L. Miao, T. Pillsbury, H. Yi, S. Kempinger, J. Hu, C. A. Heikes, P. Quarterman, et al., Phys. Rev. Res. 1, 012011 (2019).
  * Vidal et al. (2019) R. Vidal, H. Bentmann, T. Peixoto, A. Zeugner, S. Moser, C.-H. Min, S. Schatz, K. Kißner, M. Ünzelmann, C. Fornari, et al., Phys. Rev. B 100, 121104 (2019).
  * Chen et al. (2019c) Y. Chen, L. Xu, J. Li, Y. Li, H. Wang, C. Zhang, H. Li, Y. Wu, A. Liang, C. Chen, et al., Phys. Rev. X 9, 041040 (2019c).
  * Ding et al. (2021) Y. Ding, M. Zeng, Q. Zheng, J. Zhang, D. Xu, W. Chen, C. Wang, S. Chen, Y. Xie, Y. Ding, et al., Nat. Commun. 12, 5886 (2021).
  * Mazza et al. (2022) A. R. Mazza, J. Lapano, H. M. MeyerIII, C. T. Nelson, T. Smith, Y.-Y. Pai, K. Noordhoek, B. J. Lawrie, T. R. Charlton, R. G. Moore, et al., Adv. Funct. Mater. p. 2202234 (2022).
  * Giannozzi et al. (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al., J. Phys.: Condens. Matter 21, 395502 (2009).
  * Hamann (2013) D. Hamann, Phys. Rev. B 88, 085117 (2013).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Grimme et al. (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Dudarev et al. (1998) S. Dudarev, G. Botton, S. Savrasov, C. Humphreys, and A. Sutton, Phys. Rev. B 57, 1505 (1998).
  * Liechtenstein et al. (1995) A. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995).
  * Mostofi et al. (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 185, 2309 (2014).
  * Wu et al. (2018) Q. Wu, S. Zhang, H.-F. Song, M. Troyer, and A. A. Soluyanov, Comput. Phys. Commun. 224, 405 (2018).
  * Soluyanov and Vanderbilt (2011) A. A. Soluyanov and D. Vanderbilt, Phys. Rev. B 83, 035108 (2011).
  * Bennett et al. (2022) M. C. Bennett, J. Ahn, A. Pham, G. Wang, P. Ganesh, and J. T. Krogel, arXiv preprint arXiv:2204.04268 (2022).
  * Annaberdiyev et al. (2018) A. Annaberdiyev, G. Wang, C. A. Melton, M. C. Bennett, L. Shulenburger, and L. Mitas, J. Chem. Phys. 149, 134108 (2018).
  * Metz et al. (2000) B. Metz, H. Stoll, and M. Dolg, J. Chem. Phys. 113, 2563 (2000).
  * Stoll et al. (2002) H. Stoll, B. Metz, and M. Dolg, J. Comput. Chem. 23, 767 (2002).
  * Bennett et al. (2017) M. C. Bennett, C. A. Melton, A. Annaberdiyev, G. Wang, L. Shulenburger, and L. Mitas, J. Chem. Phys. 147, 224106 (2017).
  * Umrigar et al. (2007) C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hennig, Phys. Rev. Lett. 98, 110201 (2007).
  * Kim et al. (2018) J. Kim, A. D. Baczewski, T. D. Beaudet, A. Benali, M. C. Bennett, M. A. Berrill, N. S. Blunt, E. J. L. Borda, M. Casula, D. M. Ceperley, et al., J. Phys.: Condens. Matter 30, 195901 (2018).
  * Kent et al. (2020) P. R. C. Kent, A. Annaberdiyev, A. Benali, M. C. Bennett, E. J. Landinez Borda, P. Doak, H. Hao, K. D. Jordan, J. T. Krogel, I. Kylänpää, et al., J. Chem. Phys. 152, 174105 (2020).
  * Ding et al. (2020) L. Ding, C. Hu, F. Ye, E. Feng, N. Ni, and H. Cao, Phys. Rev. B 101, 020412 (2020).
  * Blaiszik et al. (2016) B. Blaiszik, K. Chard, J. Pruyne, R. Ananthakrishnan, S. Tuecke, and I. Foster, JOM 68, 2045 (2016), ISSN 1543-1851.
  * Blaiszik et al. (2019) B. Blaiszik, L. Ward, M. Schwarting, J. Gaff, R. Chard, D. Pike, K. Chard, and I. Foster, MRS Communications 9, 1125 (2019), ISSN 2159-6867.


YITP-23-121

A new perspective on thermal transition in QCD

Masanori Hanadaa, Hiroki Ohata,b Hidehiko Shimada,b and Hiromasa Watanabeb

aSchool of Mathematical Sciences, Queen Mary University of London
Mile End Road, London, E1 4NS, United Kingdom

bYukawa Institute for Theoretical Physics, Kyoto University
Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan

Abstract

Motivated by the picture of partial deconfinement developed in recent years
for large-$N$ gauge theories, we propose a new way of analyzing and
understanding thermal phase transition in QCD. We find nontrivial support for
our proposal by analyzing the WHOT-QCD collaboration’s lattice configurations
for SU(3) QCD in $3+1$ spacetime dimensions with up, down, and strange quarks.

We find that the Polyakov line (the holonomy matrix around a thermal time
circle) is governed by the Haar-random distribution at low temperatures. The
deviation from the Haar-random distribution at higher temperatures can be
measured via the character expansion, or equivalently, via the expectation
values of the Polyakov loop defined by the various nontrivial representations
of SU(3).

We find that the Polyakov loop corresponding to the fundamental representation
and loops in the higher representation condense at different temperatures.
This suggests that there are (at least) three phases, one intermediate phase
existing in between the completely-confined and the completely-deconfined
phases. Our identification of the intermediate phase is supported also by the
condensation of instantons: by studying the instanton numbers of the WHOT-QCD
configurations, we find that the instanton condensation occurs for temperature
regimes corresponding to what we identify as the completely-confined and
intermediate phases, whereas the instantons do not condense in the completely-
deconfined phase.

Our characterization of confinement based on the Haar-randomness explains why
the Polyakov loop is a good observable to distinguish the confinement and the
deconfinement phases in QCD despite the absence of the $\mathbb{Z}_{3}$ center
symmetry.

##  1 Introduction

Confinement/deconfinement transition in gauge theory [1, 2] is an important
phenomenon that has applications ranging from the quark-gluon plasma phase in
QCD under extreme conditions to the description of black holes via
gauge/gravity duality. However, the definition of confinement and
deconfinement has been somewhat unclear in real-world QCD with three colors.
The biggest issue was the lack of the strict notion of symmetry characterizing
confinement and deconfinement. Specifically, the $\mathbb{Z}_{3}$ center
symmetry that provides us with a good characterization for pure Yang-Mills
theory (i.e., unbroken and broken center symmetry correspond to confinement
and deconfinement, respectively) does not exist for QCD due to the presence of
quarks in the fundamental representation. The Polyakov loop111We will call the
holonomy matrix (in the fundamental representation) along a thermal circle as
the Polyakov line. Its trace is the standard Polyakov loop. We will also
consider the trace of the holonomy in representations other than the
fundamental representations. We will also call them the Polyakov loops in
these representations. no longer plays the role of the order parameter
associated with the center symmetry. Nonetheless, the Polyakov loop has been
empirically used as an “order parameter” of the deconfinement transitions.

Investigations of large-$N$ gauge theories have been successful in elucidating
the nontrivial nature of QCD such as the occurrence of deconfinement
transition. It has been known that the phase transition can be detected by the
Polyakov line (i.e., the holonomy matrix whose trace gives the Polyakov loop)
in an independent way to center symmetry. 222This was analytically shown for
weakly-coupled theories on a three-sphere [3, 4, 5]. Even for theories with
center symmetry, a phase transition sits between two center-broken phases. The
reason that the Polyakov line captures the phase transition turned out to be
its connection to gauge symmetry [6].

In this letter and a companion paper [7], we explore whether and how this idea
in large-$N$ theory can be applied to the real-world QCD, i.e., SU($N=3$) QCD
with dynamical quarks. Note that, although we are considering finite $N$, the
phase transition is not prohibited because the thermodynamic limit is realized
at infinite volume.

We take a bottom-up approach in this letter, by looking into lattice data
generated by WHOT-QCD collaboration 333 WHOT-QCD collaboration [8] studies the
finite temperature QCD using the Wilson fermion. Although the mass parameters
of fermions of WHOT-QCD are not small enough to reproduce the correct meson
spectrum, they are small enough so that there is no first-order thermal phase
transition, which is expected to be the behavior of QCD at physical quark
mass[9, 10].  in Sec. 2 and then giving the interpretation based on the
connection to the theoretical understanding of large-$N$ theory in Sec. 3. The
presentation in the companion paper [7] follows a top-down approach, i.e., we
start with the large-$N$ theories, make conjectures on SU(3) QCD, and then
confirm these conjectures based on lattice QCD data. In fact, we took such a
top-down approach in our investigation.

##  2 Looking into lattice data

We consider thermal circles defined at each spatial point $\vec{x}$ and the
Polyakov lines $P_{\vec{x}}\in{\rm SU}(3)$. A technical but crucial idea which
facilitates our analysis is to consider an ensemble (probability distribution)
of the Polyakov line, counting $P_{\vec{x}}$ for each $\vec{x}$ as a sample.
In Table 1, we show a short profile of the WHOT-QCD lattice configurations at
seven different values of temperature which we use, together with our
identification of the phases associated with the temperatures which will be
explained below. The spatial volume is $32^{3}$, and hence, we obtain
$32^{3}=32768$ $P_{\vec{x}}$’s from each lattice QCD configuration.

Lattice size | Temperature | Phase
---|---|---
$4\times 32^{3}$ | 697 MeV | CD
$6\times 32^{3}$ | 464 MeV | CD
$8\times 32^{3}$ | 348 MeV |  PD or CD
$10\times 32^{3}$ | 279 MeV | PD
$12\times 32^{3}$ | 232 MeV | PD
$14\times 32^{3}$ | 199 MeV | PD
$16\times 32^{3}$ | 174 MeV |  CC or PD
Table 1: WHOT-QCD configurations of $N_{\rm f}=2+1$ QCD [8]. Lattice spacing
$a\simeq 0.07$ fm, with up and down quarks heavier than physical mass (i.e.,
$\frac{m_{\pi}}{m_{\rho}}\simeq 0.63$ is larger than the real-world value
$\frac{m_{\pi}}{m_{\rho}}\simeq 0.18$) and approximately physical strange-
quark mass. The third column is our estimate of the phase discussed in the
main text. CD, PD, and CC denote Complete Deconfinement, Partial
Deconfinement, and Complete Confinement.

At each spatial point $\vec{x}$, $P_{\vec{x}}$ is a $3\times 3$ matrix with
eigenvalues $e^{i\theta_{1}},e^{i\theta_{2}}$ and $e^{i\theta_{3}}$, with
$\theta_{1}+\theta_{2}+\theta_{3}\equiv 0$ mod $2\pi$. There are $3\times
32^{3}=98304$ eigenvalues per configuration. We can estimate the distribution
$\rho(\theta)$ ($-\pi<\theta\leq\pi$) by combining many configurations. The
results are shown in Fig 1. We compare this distribution against that would
arise from the Haar-random distribution on $\rm{SU}(3)$,

$\displaystyle\rho_{\rm
Haar}(\theta)=\frac{1}{2\pi}\left(1+\frac{2}{3}\cos(3\theta)\right).$ (1)

(For completeness, we give the derivation of this formula in appendix B.) The
plots show that $\rho(\theta)$ deviates from the Haar-random distribution for
higher temperatures $T=348~{}{\rm MeV}$, $464~{}{\rm MeV}$, but appears
indistinguishable for lower temperatures, with our naked eyes, from $\rho_{\rm
Haar}(\theta)$ at $T=174~{}{\rm MeV},232~{}{\rm MeV}$. We will discuss a
method to measure the deviation from the Haar-random distribution
systematically shortly below, which shows that the deviation decreases rapidly
toward $T=174~{}{\rm MeV}$.

The agreement with the Haar-random distribution at low temperatures is a
crucial feature that has been theoretically understood in the large-$N$
theories based on analysis focused on the redundancy of the states under gauge
transformations. 444It might not come as a surprise that the Polyakov line
obeys the Haar-randomness at strictly zero temperature; at $T=0$ the Polyakov
line is a product of infinitely many link variables, and because of this the
probability distribution of the Polyakov line may converge into the Haar-
random distribution. What is quite remarkable is that the Polyakov line obeys
the Haar-random distribution to good accuracy even at rather high temperatures
right below the deconfinement transition. We provide a short summary of this
point in Appendix A.

We can measure the deviation from the Haar-randomness quantitatively by using
character expansion. 555 We note that Polyakov, in the pioneering work [1],
pointed out the usefulness of character expansion to characterize the
deconfinement transition. Apparently, this prescient remark was not followed
up seriously before our work.  (We collect some properties of character
expansion together with explicit formulae for characters of $\rm{SU}(3)$ in
Appendix C.) We write the probability distribution of the Polyakov line
$P\in\rm{SU}(3)$ by $\rho(P)$. Let $\chi_{r}(P)$ be the character associated
with an irreducible representation $r$. By the completeness of the characters,
one can expand $\rho(P)$ in terms of $\chi_{r}(P)$ as
$\rho(P)=\sum_{r}\rho_{r}\chi_{r}(P)$, where the expansion coefficients are
$\rho_{r}=\int dP\rho(P)\chi_{r}^{*}(P)$ because of the orthonormality of the
characters $\int dP\chi_{r}(P)\chi_{r^{\prime}}^{*}(P)=\delta_{rr^{\prime}}$.
By construction, $\rho_{r}$ coincides with the expectation value of the
Polyakov loop in the representation $r$.

For the exact Haar-random distribution, $\rho(P)$ is completely dominated by
the trivial representation, i.e., $\rho_{r}$ vanishes for any nontrivial
representation $r$. Hence, the Polyakov loops in nontrivial representations
give good measures of the deviation of the Polyakov line from the Haar-random
distribution. Note that $\rho_{r}$ contains all statistical information of the
distribution of the Polyakov line, including the correlations between three
eigenvalues of the Polyakov line such as the level repulsion.

Fig. 2 and Fig. 3 show the expectation values of the Polyakov loops in several
nontrivial representations. 666Note that the expectation values of Polyakov
loops are real in the absence of the chemical potential, since
$\rho(P)=\rho(P^{\dagger})$.  These plots show, firstly, that the Polyakov
loops do disappear at $T\lesssim 174$ MeV. In particular, Fig. 3 shows that
the Polyakov loop in the fundamental representation is suppressed
exponentially in the low-temperature regime. 777We cannot exclude the
possibility of exponentially small deviation from the Haar-randomness at
nonzero temperatures which is consistent with zero with our numerical
precision. This is indeed the case for large-$N$ QCD on a small three-sphere
[11]. See also the discussion section. Secondly, the expectation values of the
Polyakov loops in higher representations become nonzero at different
temperatures, $T\gtrsim 348$ MeV.

The simplest possibility consistent with our observations is that there are
three phases: (i) $T<T_{1}\sim 174$ MeV where the Polyakov line is governed by
the Haar-random distribution, (ii) $T_{1}<T<T_{2}\sim 348\rm{MeV}$ where the
fundamental Polyakov loop is non-zero but the Polyakov loops associated with
higher representations vanish, and (iii) $T_{2}<T$ where Polyakov loop in all
representations are non-zero.

It is natural to interpret the nonzero values of the Polyakov loops in
nontrivial representations as indicating that the degrees of freedom
associated with the representations are deconfined. As such, we shall call the
regimes, (i) the completely-confined, (ii) partially-deconfined (or
equivalently, partially-confined), and (iii) completely-deconfined phases,
respectively. The identified phases are shown in Table 1.

Remarkably, further support for this identification of the phases is obtained
by studying the condensation of instantons. Namely, we find that the instanton
condensation occurs for temperature regimes corresponding to what we identify
as the completely-confined and partially-confined phases, whereas the
instanton does not condense in the completely-deconfined phase. The quantity
we use to detect the instanton condensation is the topological charge of each
lattice configuration $Q$ computed by the WHOT-QCD collaboration [12]. Because
the topological charge is sensitive to the ultraviolet cutoff, it should be
evaluated from a smeared lattice configuration, for example, by using the
gradient flow [13]. After smearing, each configuration returns an integer
value, or more precisely speaking, the histogram peaks at integer values. Fig.
4 shows the distributions of the topological charge at various temperatures.
At $T\leq 279$ MeV, we clearly see multiple peaks including the ones at $Q\neq
0$ that signal the instanton condensation. (That the peaks gradually become
lower and eventually disappear can be understood as the finite-volume effect.)
The condensation melts as the temperature goes up. At $T=348$ MeV, we observe
one sharp peak at $Q=0$ and almost vanishing peaks at $Q=\pm 1$. The
instantons cease to condense around this temperature.

Figure 1: Distribution of Polyakov line phases, $N_{t}\times 32^{3}$ lattice,
obtained from the WHOT-QCD configurations, is shown in orange. The lattice
size is $N_{t}\times 32^{3}$ ($N_{t}=6,8,12,16$, and correspondingly $T=$ 464,
348, 232, and 174 MeV.) and 599 configurations were used. These histograms are
drawn with 992 bins. Although the agreement with the Haar-random distribution
shown in the dash-dot orange line (see eq. (1)) seems to be good at $T=232$
MeV, more careful investigation reveals a small deviation and hence the onset
of partial deconfinement; see Figs. 2 and 3 and main text. At 174 MeV,
agreement with the Haar-random distribution is much better.

Figure 2:  The expectation values of characters vs. temperature for the
fundamental, adjoint, rank-2 symmetric, and rank-3 symmetric representations,
obtained from WHOT-QCD configurations. Note that the expectation values are
real.

Figure 3: The log plot of the expectation value of the Polyakov loop in
fundamental representation. It is very close to zero at the lowest temperature
in the set of configurations (174 MeV).

Figure 4: Histogram of the topological charge provided by the authors of Ref.
[12].

##  3 Comparison to the large-$N$ partial deconfinement

The observations in the previous section fit nicely within the framework of
partial deconfinement [14, 15, 16, 17] in large-$N$ theories. The close
connection between the large-$N$ partial deconfinement and the behavior of QCD
justifies the use of the terms we used to designate the phases: completely-
confined, partially-deconfined, and completely-deconfined.

To understand the meaning of partial deconfinement and Polyakov loop, the
amount of gauge redundancy plays an important role [6]. This is explained in
Appendix A. See also the companion paper [7] which explains more details.

Several researchers pointed out that the confinement/deconfinement transition
consists of two phase transitions based on the weak coupling analysis [4, 3].
A more explicit understanding of the physical interpretation and the mechanism
of the emergence of the intermediate phase was developed in a series of papers
[16, 17, 11, 6, 7]. Specifically, this phase was identified as the coexistence
of confined and deconfined degrees of freedom in the space of colors (internal
space) rather than in the usual coordinate space.

We can summarize the connection between the large-$N$ partial deconfinement
and the new perspectives discussed above as follows.

  1. 1.

The completely-confined phase is governed by the Haar-random distribution of
the Polyakov lines [6]. This is common to both large-$N$ theories and
finite-$N$ theories including QCD.

  2. 2.

The transition from the partially-deconfined phase to the completely-
deconfined phase is identified with the Gross-Witten-Wadia (GWW) transition in
the large-$N$ theory. The GWW transition can be captured by using the
expectation values of the multiply-wound Polyakov loops $u_{n}=\langle\tr
P^{n}\rangle$. In particular, after the GWW transition, all $u_{n}$ become
non-zero. 888More generally, $u_{n}$ with large $n$ may take nonzero but
exponentially suppressed values below the GWW point, as is the case for
large-$N$ QCD on a small three-sphere [5, 11]. For the QCD, we find that the
transition from the partially-deconfined phase to the completely-deconfined
phase is associated with the onsets of the Polyakov loops in the higher
representations. The Polyakov loops in the higher representations and
$u_{n}$’s (with $n\geq 1$) play analogous roles and some of them are directly
related. In Appendix C, we give examples of direct relations between $u_{n}$
and Polyakov loops in higher representations.

  3. 3.

As discussed and observed in Refs. [11, 18], it is natural to expect that the
chiral symmetry breaking takes place at the GWW point when quarks are
massless. (Intuitively, quarks in the confined sector should form a chiral
condensate; otherwise, the ’t Hooft anomaly would not be preserved.) Since the
instantons are intimately connected with the chiral symmetry, it is natural to
expect the behavior of instantons to change across the GWW transition. 999 It
was shown in Ref. [19] that the GWW transition of 2D lattice Yang-Mills theory
can be directly understood as arising from the change of the saddle points
contributing to the path integral including the contribution of the
instantons.  Since in QCD the quark mass breaks the chiral symmetry
explicitly, the instanton condensation is a natural probe for the finite-$N$
analog of the GWW transition. Indeed, we find that the phase structure
suggested by instantons is consistent with that obtained from the Polyakov
loops.

##  4 Discussion

Partial deconfinement is a generic property of gauge theories insensitive to
the details of the theories such as gauge group or matter content, at least at
large $N$. Therefore, it is natural to expect other finite-$N$ theories to
exhibit thermal phases that can be characterized in a manner similar to what
we have done for QCD. By studying them, we might be able to improve the
understanding of confinement and deconfinement. Practically, SU(3) QCD with
finite quark mass and pure Yang-Mills are the most tractable targets because
many sets of lattice configurations are available for SU(3) QCD and the
simulation of pure Yang-Mills is not costly.

It is clearly important to further verify the new perspective we advocated in
this letter. For example, studying the Polyakov loop in larger
representations, identifying the transition temperature more precisely, and
identifying the order of the phase transitions are important future problems.

The first ‘transition’ at $T_{1}$ may well be a crossover. For example, the
Polyakov loops may have small nonzero values that become exponentially small
when the sizes of the representations are increased. 101010In the case of
large-$N$ QCD on a small three-sphere, this is indeed the case [5, 11].  A
small deviation from the Haar-randomness could arise which we may attribute to
the contributions of hadron gas. The contributions are exponentially
suppressed for heavy particles.

It would be more natural to expect a transition with non-analyticity at
$T_{2}$, given the connection to the condensation of instantons. This can be
the case even if the Polyakov loops in large representations are not exactly
zero at $T<T_{2}$; a natural possibility would be that there is a transition
between the exponential decay at $T<T_{2}$ and the power-law decay at
$T>T_{2}$.

It is remarkable that the Polyakov line obeys Haar-randomness to good accuracy
for low temperatures. We have seen that the deviation from the Haar-randomness
can be measured via the character expansion. In this letter, we focused on the
one-point function of the Polyakov loops. It is interesting to consider
whether one can also study other observables from the point of view of the
Haar-randomness and the deviation from it.

An obstacle to generalizing partial deconfinement to finite $N$ had been the
meaning of the size of the deconfined sector $M$ ($0\leq M\leq N$). Even in
the large-$N$ limit, $M$ is not literally an integer. It could have an
uncertainty of order $N^{0}$ which is negligible at large $N$. Admittedly,
such an ambiguity makes the use of “$M$” very subtle at $N=3$. To circumvent
this issue, we avoided the use of $M$ and relied on characters (the Polyakov
loops in various representations). The use of the character also has the
advantage of being manifestly gauge invariant. Although in large $N$-theory
the use of the parameter $M$ is shown to have gauge invariant meaning, it may
be worthwhile to revisit the analysis of the partial deconfinement in the
large-$N$ theory using the character expansion. Character expansion played an
important role in the large-$N$ theory, for example in Refs. [20, 21]. A
recent work [22] employs the character expansion to study the deconfinement
transition in large-$N$ theory.

In this letter, we have explored the finite-$N$ counterpart of the large-$N$
partial deconfinement. The original motivation for partial deconfinement [14]
was to study black hole geometry via holography. Since $1/N$ corrections
should play crucial roles in black hole physics, we hope that our work would
be useful in obtaining important intuition into quantum gravitational
phenomena such as black hole evaporation from the QFT side.

Data management

No additional research data beyond the data presented and cited in this work
are needed to validate the research findings in this work.

Acknowledgement

We would like to thank the members of the WHOT-QCD collaboration, including
Shinji Ejiri, Kazuyuki Kanaya, Masakiyo Kitazawa, and Takashi Umeda, for
providing us with their lattice configurations and many plots and having
stimulating discussions with us. The analysis of topological charge in Ref.
[12] was led by Yusuke Taniguchi who passed away in 2022. Kazuyuki Kanaya
collected the data and plots created by Yusuke Taniguchi for us. We deeply
thank Yusuke Taniguchi and Kazuyuki Kanaya. We also thank Sinya Aoki, Hidenori
Fukaya, Vaibhav Gautam, Yui Hayashi, Jack Holden, Seok Kim, Tamas Kovacs,
Atsushi Nakamura, Marco Panero, Robert Pisarski, Enrico Rinaldi, Yuya
Tanizaki, and Jacobus Verbaarschot for discussions and comments.

This work was supported by the Japan Lattice Data Grid (JLDG) constructed over
the SINET5 of NII. M. H. thanks for the STFC grant ST/R003599/1. H. O. is
supported by a Grant-in-Aid for JSPS Fellows (Grant No.22KJ1662). H. S. is
supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI
Grant number 21H05182. H. W. is supported by the Japan Society for the
Promotion of Science (JSPS) KAKENHI Grant number 22H01218.

##  Appendix A Polyakov loop and the amount of redundancy at $N=\infty$

The purpose of this appendix is to give a short summary of the essential
mechanism of the large-$N$ partial deconfinement. In particular, we explain
the relevance of the redundancy of states under gauge transformations and its
relation to the Polyakov line.

That the amount of gauge redundancy has important consequences is not really
new. In fact, it has been known for a century since the theoretical discovery
of Bose-Einstein condensation [23], although the connection to Polyakov line
and confinement was pointed out only recently [6]. To see how the Polyakov
line and gauge redundancy are related, we consider three theories with
increasing levels of complexity: $N$ indistinguishable bosons, SU($N$)
Hermitian matrix model, and SU($N$) QCD.

####  $N$ indistinguishable bosons and Bose-Einstein condensation

To describe $N$ indistinguishable bosons in $\mathbb{R}^{3}$, we use
coordinate operators $\hat{\vec{x}}_{i}=(\hat{x}_{i},\hat{y}_{i},\hat{z}_{i})$
and momentum operators
$\hat{\vec{p}}_{i}=(\hat{p}_{x,i},\hat{p}_{y,i},\hat{p}_{z,i})$, where
$i=1,2,\cdots,N$ labels bosons. The Hamiltonian $\hat{H}$ is invariant under
the permutation of the labels, e.g.,
$\hat{H}=\sum_{i=1}^{N}(\frac{1}{2m}\hat{\vec{p}}_{i}^{2}+\frac{m\omega^{2}}{2}\hat{\vec{x}}_{i}^{2})$
in the weak-coupling limit. That the bosons are indistinguishable means the SN
permutation symmetry is gauged.

We can use the coordinate eigenstates
$\ket{x}=\ket{\vec{x}_{1},\cdots,\vec{x}_{N}}$ that satisfy
$\hat{\vec{x}}_{i}\ket{x}=\vec{x}_{i}\ket{x}$ to describe quantum states. The
coordinate eigenstates span the extended Hilbert space ${\cal H}_{\rm ext}$
that contains SN-non-singlets:

$\displaystyle\mathcal{H}_{\rm
ext}=\textrm{Span}\\{\ket{x}|x\in\mathbb{R}^{3N}\\}.$ (2)

For a permutation $\sigma\in\text{S}_{N}$, we define $\hat{\sigma}$ by
$\hat{\sigma}\ket{\vec{x}_{1},\cdots,\vec{x}_{N}}=\ket{\vec{x}_{\sigma(1)},\cdots,\vec{x}_{\sigma(N)}}$.
From this, we can define the projection operator
$\hat{\pi}=\frac{1}{N!}\sum_{\sigma\in\text{S}_{N}}\hat{\sigma}$ that maps
$\mathcal{H}_{\rm ext}$ to the SN-invariant Hilbert space $\mathcal{H}_{\rm
inv}$. Canonical partition function at temperature $T$ can be written in two
ways, in terms of $\mathcal{H}_{\rm ext}$ and $\mathcal{H}_{\rm inv}$:

$\displaystyle Z(T)=\mathrm{Tr}_{\mathcal{H}_{\rm
inv}}(e^{-\hat{H}/T})=\mathrm{Tr}_{\mathcal{H}_{\rm
ext}}(\hat{\pi}e^{-\hat{H}/T})=\frac{1}{N!}\sum_{\sigma\in\mathrm{S}_{N}}\mathrm{Tr}_{\mathcal{H}_{\rm
ext}}(\hat{\sigma}e^{-\hat{H}/T})\,.$ (3)

To see how states in $\mathcal{H}_{\rm inv}$ and $\mathcal{H}_{\rm ext}$ are
related, let us consider the weak-coupling limit and consider a product of
one-particle states, $\ket{\Phi}$, whose wave function is written as
$\bra{x}\ket{\Phi}=\prod_{i=1}^{N}\phi_{i}(\vec{x}_{i})$. The symmetric group
SN acts on $\ket{\Phi}$ as
$\bra{x}\hat{\sigma}\ket{\Phi}=\prod_{i=1}^{N}\phi_{i}(\vec{x}_{\sigma(i)})$.
If all one-particle states are the same, i.e., $\phi_{1}=\cdots=\phi_{N}$, the
product state is invariant under SN. On the other hand, if
$\phi_{1},\cdots,\phi_{M}$ are all different but $\phi_{M+1}=\cdots=\phi_{N}$,
then the product state is invariant under
$\mathrm{S}_{N-M}\subset\mathrm{S}_{N}$. Such an unbroken symmetry acting on
$\mathcal{H}_{\rm ext}$ leads to an enhancement factor in the partition
function. Specifically, when the overlap of different one-particle states can
be neglected, we obtain

$\displaystyle\sum_{\sigma\in\mathrm{S}_{N}}\bra{\Phi}\hat{\sigma}\ket{\Phi}=(N-M)!\,.$
(4)

This enhancement factor is responsible for the Bose-Einstein condensation.
Many particles fall into the one-particle ground state assisted by this
enhancement factor. From this, one can learn that less redundant states are
preferred at low energy.

We will next discuss the SU($N$) Hermitian matrix model and will see that the
permutation $\sigma$ is nothing but the Polyakov line. This implies that the
typical Polyakov line in the path integral is determined in such a way that it
leaves typical states dominating partition function invariant, i.e.,
$\hat{\sigma}\ket{\rm typical}=\ket{\rm typical}$. We will elaborate on this
statement shortly.

#### SU($N$) Hermitian matrix model and color confinement

Let us consider a bosonic SU($N$) gauged $D$-matrix model ($D\geq 2$) with the
Hamiltonian

$\displaystyle\hat{H}={\rm
tr}\left(\frac{1}{2}\sum_{I=1}^{D}\hat{P}_{I}^{2}+V(\hat{X})\right)\,,$ (5)

where ${\rm tr}$ means the trace as $N\times N$ matrix and $V(\hat{X})$ is a
potential term such as
$V(\hat{X})=-\frac{g^{2}}{4}[\hat{X}_{I},\hat{X}_{J}]^{2}$. Each $\hat{X}_{I}$
has $N^{2}$ components $\hat{X}_{I,ij}$, where $i,j=1,2,\cdots,N$, that
satisfy the Hermiticity condition $(\hat{X}_{I,ij})^{\dagger}=\hat{X}_{I,ji}$.
We do not impose the traceless condition. The operator $\hat{P}_{I,ij}$ is the
conjugate momentum of $\hat{X}_{I,ji}$. They satisfy the canonical commutation
relation $[\hat{X}_{I,ij},\hat{P}_{J,kl}]=i\delta_{IJ}\delta_{il}\delta_{jk}$.

The Hamiltonian is invariant under the adjoint action of SU($N$) defined by

$\displaystyle\hat{X}_{I,ij}\to(U\hat{X}_{I}U^{-1})_{ij}=\sum_{k,l}U_{ij}\hat{X}_{I,kl}U^{-1}_{lj},\qquad\hat{P}_{I,ij}\to(U\hat{P}_{I}U^{-1})_{ij}=\sum_{k,l}U_{ik}\hat{P}_{I,kl}U^{-1}_{lj}.$
(6)

We can use the extended Hilbert space with SU($N$) non-singlet states, ${\cal
H}_{\rm ext}$. The extended space is spanned by the coordinate eigenstates
$\ket{X}$ that satisfy $\hat{X}_{I,ij}\ket{X}=X_{I,ij}\ket{X}$:

$\displaystyle\mathcal{H}_{\rm
ext}=\textrm{Span}\\{\ket{X}|X\in\mathbb{R}^{DN^{2}}\\}.$ (7)

Note that $X$ consists of $DN^{2}$ real numbers
$X_{I=1,\cdots,D}^{\alpha=1,\cdots,N^{2}}$. Gauge transformation acts on
$\mathcal{H}_{\rm ext}$ as $\ket{X}\to\ket{U^{-1}XU}$. By using the
SU($N$)-invariant Hilbert space ${\cal H}_{\rm inv}$, the canonical partition
function at temperature $T$ can be written as [6]

$\displaystyle Z(T)={\rm Tr}_{\mathcal{H}_{\rm
inv}}\left(e^{-\hat{H}/T}\right)=\mathrm{Tr}_{\mathcal{H}_{\rm
ext}}(\hat{\pi}e^{-\hat{H}/T})=\frac{1}{\textrm{Vol}G}\int_{G}dg\mathrm{Tr}_{\mathcal{H}_{\rm
ext}}(\hat{g}e^{-\hat{H}/T})\,,$ (8)

where $G=\mathrm{SU}(N)$ and
$\hat{\pi}\equiv\frac{1}{\textrm{vol}G}\int_{G}dg\,\hat{g}$ is a projection
operator from ${\cal H}_{\rm ext}$ to ${\cal H}_{\rm inv}$.

Figure 5: Partial deconfinement in the matrix model. The $M\times M$-block
shown in red is deconfined. This figure is taken from Ref. [17].

The standard technique to rewrite the Hamiltonian formulation to the path-
integral formulation tells us that $g\in G$ is the Polyakov line [6]. By
comparing (8) with (3), we see that $\sigma\in\mathrm{S}_{N}$ corresponds to
the Polyakov line. The counterpart of the Bose-Einstein condensation is
confinement [6]. In general, SU($M$) subgroup of SU($N$) can be deconfined
[14, 15, 16, 17], leaving $\mathrm{SU}(N-M)\subset\mathrm{SU}(N)$ as a
symmetry in $\mathcal{H}_{\rm ext}$ that leads to the enhancement factor
$\mathrm{Vol}(\mathrm{SU}(N-M))\sim e^{(N-M)^{2}}$ [6]. The value of $M$
depends on the energy $E$ in a nontrivial manner.

#### Distribution of Polyakov line phases in the large-$N$ limit

Let us further focus on typical states dominating thermodynamics that are
invariant under $\textrm{SU}(N-M)\subset\textrm{SU}(N)$. We can fix the gauge
in such a way that the SU($M$)-deconfined sector is the upper-left $M\times
M$-block as in Fig. 5. This choice of embedding of SU($M$) into SU($N$) fixes
SU($N$) down to $\mathrm{SU}(M)\times\mathrm{SU}(N-M)$. Then, the Polyakov
line $\mathcal{P}$ takes the following form:

$\displaystyle\mathcal{P}=\left(\begin{array}[]{cc}\mathcal{P}_{\rm dec}&0\\\
0&\mathcal{P}_{\rm con}\end{array}\right)\,.$ (11)

The eigenvalues of $\mathcal{P}_{\rm dec}$ and $\mathcal{P}_{\rm con}$ give
the phases $\theta_{1},\cdots,\theta_{M}$ and
$\theta_{M+1},\cdots,\theta_{N}$, respectively. From them, we can determine
the distribution of the phases $\rho_{\rm dec}(\theta)$ and $\rho_{\rm
con}(\theta)$. The latter is constant,

$\displaystyle\rho_{\rm con}(\theta)=\frac{1}{2\pi}\,.$ (12)

This is because $\mathcal{P}_{\rm con}$ can be any element of SU($N-M$)
(specifically, $\mathcal{P}_{\rm con}$ dominating path integral is distributed
following the Haar measure) and the generic phase distribution in this part is
uniform in the limit of $N-M\to\infty$. The former is not uniform and its
smallest value is zero. The full distribution is [6, 16, 17]

$\displaystyle\rho(\theta)=\frac{1}{2\pi}\cdot\left(1-\frac{M}{N}\right)+\frac{M}{N}\cdot\rho_{\rm
dec}(\theta)\,.$ (13)

Therefore, the constant offset provides us with a way to measure the amount of
gauge redundancy:

$\displaystyle\mathrm{\ Constant\
offset}=\frac{1}{2\pi}\left(1-\frac{M}{N}\right)\,.$ (14)

This relation holds for the system of indistinguishable bosons as well.

#### SU($N$) QCD

Essentially the same mechanism works for large-$N$ QCD. Although we leave the
detail to the companion paper [7], the essence is easy to understand: If the
fields are slowly varying up to gauge transformations (which must be the case
at low energy), in a small enough spatial volume we can reduce the theory to a
matrix model (i.e., a theory without spatial extent), and then, the same
mechanism analogous to BEC works.

##  Appendix B SU($N$) Haar-random distribution

In this appendix, we show the analytic formula for the Haar-random
distribution of the phases in SU($N$),

$\displaystyle\rho_{\rm
Haar}(\theta)=\frac{1}{2\pi}\left(1-(-1)^{N}\cdot\frac{2}{N}\cos(N\theta)\right)\,,$
(15)

for $N=3$. See Ref. [7] for other values of $N$.

We start with

$\quantity[\differential
U]=\rho(\theta_{1},\cdots,\theta_{N})\prod_{j=1}^{N}\differential\theta_{j}\,,$
(16)

where

$\displaystyle\rho(\theta_{1},\cdots,\theta_{N})=C\prod_{j<k}\sin^{2}\left(\frac{\theta_{j}-\theta_{k}}{2}\right)\times\sum_{n=-\infty}^{\infty}\delta\left(\sum_{j=1}^{N}\theta_{j}-2\pi
n\right)\,.$ (17)

Here, $C$ is the normalization factor and
$\prod_{j<k}\sin^{2}\left(\frac{\theta_{j}-\theta_{k}}{2}\right)$ is the
Vandermonde determinant. Our task is to integrate out
$\theta_{2},\cdots,\theta_{N}$ and estimate

$\displaystyle\rho(\theta_{1})=\int\differential\theta_{2}\cdots\int\differential\theta_{N}\
\rho(\theta_{1},\theta_{2},\cdots,\theta_{N})\,.$ (18)

For $N=3$, we use $\theta_{3}=-(\theta_{1}+\theta_{2})$ up to the shift by
$2\pi n$. By substituting this into (17), we obtain

$\displaystyle\rho(\theta_{1},\theta_{2})$
$\displaystyle=C\sin^{2}\left(\frac{\theta_{1}-\theta_{2}}{2}\right)\sin^{2}\left(\frac{\theta_{1}+2\theta_{2}}{2}\right)\sin^{2}\left(\frac{\theta_{2}+2\theta_{1}}{2}\right)$
$\displaystyle=\frac{C}{8}\left(1-\cos(\theta_{1}-\theta_{2})\right)\left(1-\cos(\theta_{1}+2\theta_{2})\right)\left(1-\cos(2\theta_{1}+\theta_{2})\right)$
$\displaystyle=\frac{C}{8}\left(\frac{3}{4}+\frac{1}{2}\cos(3\theta_{1})+\cdots\right)\,,$
(19)

where $\cdots$ are terms like $\cos(\theta_{1}-\theta_{2})$ or
$\cos(\theta_{1}+2\theta_{2})$ which disappear after $\theta_{2}$ is
integrated. By integrating $\theta_{2}$, we obtain

$\rho(\theta_{1})=\int_{-\pi}^{\pi}\differential\theta_{2}\,\rho(\theta_{1},\theta_{2})\propto
1+\frac{2}{3}\cos(3\theta_{1})\,.$ (20)

We can fix the overall factor from
$\int_{-\pi}^{\pi}\differential\theta_{1}\,\rho(\theta_{1})=1$.

##  Appendix C Characters of SU(3) group

The character of $G$ is defined by the trace of the representation matrix
$R_{r}$ in the representation $r$ as

$\chi_{r}(g):=\tr R_{r}(g).$ (21)

For the irreducible representations, the orthonormal condition,

$\frac{1}{\mathrm{Vol}(G)}\int_{G}\mathrm{d}g\,\chi_{r}(g)\left(\chi_{s}(g)\right)^{\ast}=\delta_{rs},$
(22)

is satisfied.

When we consider $G=$ SU(3), we can identify $g$ and its representation in the
fundamental representation: $g$ is $3\times 3$ unitary matrix with a
determinant equal to 1. In our context $g$ is identified with the Polyakov
line $P$. The SU(3) characters are functions only with respect to the
eigenvalues of $g$, denoted by $\lambda_{j}$ $(j=1,2,3)$. Namely,

$\chi_{r}(g)=\chi_{r}(\\{\lambda\\}).$ (23)

Note that $\lambda_{j}=e^{i\theta_{j}}$ in terms of the Polyakov line phase
$\theta_{j}$, and $\lambda_{1}\lambda_{2}\lambda_{3}=1$ due to the condition
$\det g=1$.

###  C.1 List of characters for irreducible representations

The character is a symmetric polynomial in $\lambda$’s. As is well-known,
irreducible representations of $\rm{SU}(3)$ can be labelled by Young diagram
with two rows. In the following, the notation $(n,m)$ represents a Young
tableau that has $n$ and $m$ boxes in the first and second rows, respectively.

(0,0)=1: trivial

$\chi_{\rm trivial}\equiv\chi_{(0,0)}=1\,.$ (24)

(1,0)=3: fundamental

$\chi_{\rm fund.}\equiv\chi_{(1,0)}=\lambda_{1}+\lambda_{2}+\lambda_{3}\,.$
(25)

(1,1)=$\bar{\bf 3}$: rank-two anti-symmetric, or anti-fundamental

$\chi_{\rm
antisym.}\equiv\chi_{(1,1)}=\lambda_{1}\lambda_{2}+\lambda_{2}\lambda_{3}+\lambda_{3}\lambda_{1}=\lambda_{1}^{-1}+\lambda_{2}^{-1}+\lambda_{3}^{-1}=(\chi_{\rm
fund.})^{\ast}\,.$ (26)

(2,0)=6: rank-two symmetric

$\chi_{(2,0)}=\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}+\lambda_{1}\lambda_{2}+\lambda_{2}\lambda_{3}+\lambda_{3}\lambda_{1}\,.$
(27)

(2,1)=8: adjoint

$\chi_{\rm adj.}\equiv\chi_{(2,1)}=2+\sum_{a\neq
b}\lambda_{a}\lambda_{b}^{-1}\,.$ (28)

(3,0)=10: rank-three symmetric

$\displaystyle\chi_{\rm 3\mathchar 45\relax
sym.}\equiv\chi_{(3,0)}=1+\lambda_{1}^{3}+\lambda_{2}^{3}+\lambda_{3}^{3}+\lambda_{1}^{2}\lambda_{2}+\lambda_{1}^{2}\lambda_{3}+\lambda_{2}^{2}\lambda_{1}+\lambda_{2}^{2}\lambda_{3}+\lambda_{3}^{2}\lambda_{1}+\lambda_{3}^{2}\lambda_{2}.$
(29)

###  C.2 Multiply-wound Polyakov loops and characters

The multiply-wound loops are expressed as

$\displaystyle
u_{n}\equiv\frac{1}{N}\mathrm{Tr}P^{n}=\frac{1}{N}\sum_{j=1}^{N}\lambda_{j}^{n}\,.$
(30)

For the $N=3$, the following relations hold:

$\displaystyle 3u_{1}$ $\displaystyle=\chi_{(1,0)}=\chi_{\rm fund.}\,,$ (31)
$\displaystyle 3u_{2}$ $\displaystyle=\chi_{(2,0)}-\chi_{(1,1)}=\chi_{\rm
2\mathchar 45\relax sym.}-(\chi_{\rm fund.})^{\ast}\,,$ (32) $\displaystyle
3u_{3}$ $\displaystyle=\chi_{(3,0)}-\chi_{(2,1)}+\chi_{(0,0)}=\chi_{\rm
3\mathchar 45\relax sym.}-\chi_{\rm adj.}+\chi_{\rm trivial}\,,$ (33)
$\displaystyle 3u_{4}$
$\displaystyle=\chi_{(4,0)}-\chi_{(3,1)}+\chi_{(1,0)}\,,$ (34) $\displaystyle
3u_{5}$ $\displaystyle=\chi_{(5,0)}-\chi_{(4,1)}+\chi_{(2,0)}\,,$ (35)
$\displaystyle 3u_{6}$
$\displaystyle=\chi_{(6,0)}-\chi_{(5,1)}+\chi_{(3,0)}\,,$ (36)
$\displaystyle\qquad\qquad\quad\vdots$

Note that only $u_{3}$ contains $\chi_{\mathrm{trivial}}$.

Nonzero expectation values of $\chi_{r}$ for nontrivial representations $r$
characterize the deviation from the SU(3) Haar-random distribution. We
emphasize that the Haar-randomness provides stronger restriction than that by
the $\mathbb{Z}_{3}$ center symmetry: center symmetry allows nonzero values of
$u_{3n}$. The Haar-randomness, however, leads to $\chi_{\mathrm{trivial}}=1$
and $\chi_{r}=0$ for any other representation $r$, which allows only $u_{3}$
to be nonzero.

## References

  * [1] A. M. Polyakov, “Thermal Properties of Gauge Fields and Quark Liberation,” Phys. Lett. B 72 (1978) 477–480.
  * [2] L. Susskind, “Lattice Models of Quark Confinement at High Temperature,” Phys. Rev. D 20 (1979) 2610–2618.
  * [3] B. Sundborg, “The Hagedorn transition, deconfinement and N=4 SYM theory,” Nucl. Phys. B 573 (2000) 349–363, arXiv:hep-th/9908001.
  * [4] O. Aharony, J. Marsano, S. Minwalla, K. Papadodimas, and M. Van Raamsdonk, “The Hagedorn - deconfinement phase transition in weakly coupled large N gauge theories,” Adv. Theor. Math. Phys. 8 (2004) 603–696, arXiv:hep-th/0310285.
  * [5] H. J. Schnitzer, “Confinement/deconfinement transition of large N gauge theories with N(f) fundamentals: N(f)/N finite,” Nucl. Phys. B 695 (2004) 267–282, arXiv:hep-th/0402219.
  * [6] M. Hanada, H. Shimada, and N. Wintergerst, “Color confinement and Bose-Einstein condensation,” JHEP 08 (2021) 039, arXiv:2001.10459 [hep-th].
  * [7] M. Hanada, H. Ohata, and H. Watanabe, “On thermal transition in qcd,” To appear (2023) .
  * [8] WHOT-QCD Collaboration, T. Umeda, S. Aoki, S. Ejiri, T. Hatsuda, K. Kanaya, Y. Maezawa, and H. Ohno, “Equation of state in 2+1 flavor QCD with improved Wilson quarks by the fixed scale approach,” Phys. Rev. D 85 (2012) 094508, arXiv:1202.4719 [hep-lat].
  * [9] Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabo, “The QCD transition temperature: Results with physical masses in the continuum limit,” Phys. Lett. B 643 (2006) 46–54, arXiv:hep-lat/0609068.
  * [10] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, “The Order of the quantum chromodynamics transition predicted by the standard model of particle physics,” Nature 443 (2006) 675–678, arXiv:hep-lat/0611014.
  * [11] M. Hanada and B. Robinson, “Partial-Symmetry-Breaking Phase Transitions,” Phys. Rev. D 102 no. 9, (2020) 096013, arXiv:1911.06223 [hep-th].
  * [12] Y. Taniguchi, K. Kanaya, H. Suzuki, and T. Umeda, “Topological susceptibility in finite temperature ( 2+1 )-flavor QCD using gradient flow,” Phys. Rev. D 95 no. 5, (2017) 054502, arXiv:1611.02411 [hep-lat].
  * [13] M. Lüscher, “Properties and uses of the Wilson flow in lattice QCD,” JHEP 08 (2010) 071, arXiv:1006.4518 [hep-lat]. [Erratum: JHEP 03, 092 (2014)].
  * [14] M. Hanada and J. Maltz, “A proposal of the gauge theory description of the small Schwarzschild black hole in AdS${}_{5}\times$S5,” JHEP 02 (2017) 012, arXiv:1608.03276 [hep-th].
  * [15] D. Berenstein, “Submatrix deconfinement and small black holes in AdS,” JHEP 09 (2018) 054, arXiv:1806.05729 [hep-th].
  * [16] M. Hanada, G. Ishiki, and H. Watanabe, “Partial Deconfinement,” JHEP 03 (2019) 145, arXiv:1812.05494 [hep-th]. [Erratum: JHEP 10, 029 (2019)].
  * [17] M. Hanada, A. Jevicki, C. Peng, and N. Wintergerst, “Anatomy of Deconfinement,” JHEP 12 (2019) 167, arXiv:1909.09118 [hep-th].
  * [18] M. Hanada, J. Holden, M. Knaggs, and A. O’Bannon, “Global symmetries and partial confinement,” JHEP 03 (2022) 118, arXiv:2112.11398 [hep-th].
  * [19] P. V. Buividovich, G. V. Dunne, and S. N. Valgushev, “Complex Path Integrals and Saddles in Two-Dimensional Gauge Theory,” Phys. Rev. Lett. 116 no. 13, (2016) 132001, arXiv:1512.09021 [hep-th].
  * [20] D. J. Gross and W. Taylor, “Two-dimensional QCD is a string theory,” Nucl. Phys. B 400 (1993) 181–208, arXiv:hep-th/9301068.
  * [21] D. J. Gross and W. Taylor, “Twists and Wilson loops in the string theory of two-dimensional QCD,” Nucl. Phys. B 403 (1993) 395–452, arXiv:hep-th/9303046.
  * [22] D. Berenstein and K. Yan, “The endpoint of partial deconfinement,” arXiv:2307.06122 [hep-th].
  * [23] A. Einstein, “Quantentheorie des einatomigen idealen gases,” S-B Preuss. Akad. Berlin (1924) .