# Canonical partition function and center symmetry breaking in finite density
lattice gauge theories

Shinji Ejiri  Department of Physics, Niigata University, Niigata 950-2181,
Japan

(November 28, 2022)

###### Abstract

We study the nature of the phase transition of lattice gauge theories at high
temperature and high density by focusing on the probability distribution
function, which represents the probability that a certain density will be
realized in a heat bath. The probability distribution function is obtained by
creating a canonical partition function fixing the number of particles from
the grand partition function. However, if the $Z_{3}$ center symmetry, which
is important for understanding the finite temperature phase transition of
$SU(3)$ lattice gauge theory, is maintained on a finite lattice, the
probability distribution function is always zero, except when the number of
particles is a multiple of three. For $U(1)$ lattice gauge theory, this
problem is more serious. The probability distribution becomes zero when the
particle number is nonzero. This problem is essentially the same as the
problem that the expectation value of the Polyakov loop is always zero when
calculating with finite volume. In this study, we propose a solution to this
problem. We also propose a method to avoid the sign problem, which is an
important problem at finite density, using the center symmetry. In the case of
$U(1)$ lattice gauge theory with heavy fermions, numerical simulations are
actually performed, and we demonstrate that the probability distribution
function at a finite density can be calculated by the method proposed in this
study. Furthermore, the application of this method to QCD is discussed.

##  I Introduction

Various results have been obtained so far in heavy-ion collision experiments
to generate quark-gluon plasma. In particular, recently, there has been a
great deal of interest in the phase transition of QCD in the high-density
region. One of the most interesting topics is the verification that the first-
order phase transition appears at high density, which is expected from the
phenomenological discussion. Since the critical point at which the low-density
crossover changes to the first-order phase transition is the second-order
phase transition, the fluctuation of the thermodynamic quantity of the heat
bath increases as the particle density approaches the critical point.

In order to investigate the existence of such a critical point, the
measurement of event-by-event fluctuations for conserved quantities such as
net quark number and charge amount is attracting attention. Furthermore,
higher-order cumulants such as kurtosis and skewness have characteristic
behaviors near the critical point, so they are also important physical
quantities Ejiri:2005wq ; Stephanov:2008qz ; Asakawa:2009aj ; Cheng:2008zh .
On the other hand, in the numerical calculation of lattice QCD, which is the
first-principles calculation, variance, skewness, and kurtosis can be
calculated by the Taylor expansion method in the low-density region
Allton:2003vx ; Allton:2005gk ; Ejiri:2009hq ; DElia:2016jqh ; Bazavov:2017dus
; Borsanyi:2018grb ; Bazavov:2020bjn ; Bellwied:2021nrt . These are
information about the shape of the probability distribution function
(histogram) of the quark number or the charge.

The probability distribution function of the quark number can be obtained by
computing the canonical partition function Hasenfratz:1991ax ;
Kratochvila:2005mk ; deForcrand:2006ec ; Ejiri:2008xt ; Alexandru:2005ix ;
Li:2010qf ; Li:2011ee ; Danzer:2012vw ; Fukuda:2015mva ; Nakamura:2015jra ;
Bornyakov:2016wld ; Wakayama:2018wkc . In this study, we calculate the
probability distribution function of quark number. The probability
distribution function is a function that expresses the probability that a
physical quantity will be a certain value. Then, we aim to establish a method
to determine the critical chemical potential (density) by numerical
simulations of lattice QCD. If there is a first-order phase transition, two
states appear with equal probability at the phase transition point. Therefore,
if we calculate the probability distribution function and investigate its
shape, we can determine whether the phase transition is first order, second
order, or crossover.

The numerical calculation of the canonical partition function was performed in
Ref. Ejiri:2008xt  by a saddle point approximation in two flavor lattice QCD.
In the high-density region, the probability distribution function at the phase
transition point has two peaks, indicating that the first-order phase-
transition region exists. However, since it is a finite density system, the
sign problem arises in the calculation of the canonical partition function. In
Ref. Ejiri:2008xt , the sign problem is avoided by introducing an
approximation that assumes that the probability distribution of the complex
phase of the quark determinant, which causes the sign problem, is a Gaussian
function. This approximation can be justified in the high-temperature phase,
but in the low-temperature phase, the validity of the approximation is only
that the actual numerical result is similar to the approximation function.

Moreover, the calculation of Ref. Ejiri:2008xt  is implicitly limited to the
case where the number of quarks is a multiple of three. This is because, when
calculated straightforwardly, the probability distribution function becomes
zero except for multiples of three due to the $Z_{3}$-center symmetry. The
$Z_{3}$-center symmetry is an important symmetry for understanding the
deconfinement phase transition of QCD as a spontaneous symmetry breaking. If
the symmetry is strictly maintained under the $Z_{3}$ center transformation,
only states where the net quark number is zero or a multiple of three are
realized. It is consistent with the existence of only mesons (quark and
antiquark) and baryons (three quarks) in the confined phase where symmetry is
maintained. The deconfinement phase is realized by the spontaneous symmetry
breaking, which allows the quark to exist alone. However, the symmetry is not
broken in the actual numerical calculation with finite volume. In other words,
as long as there is symmetry, the number of quarks will always be a multiple
of three.

In order to discuss phase transitions due to spontaneous symmetry breaking by
a correct procedure, it is necessary to add an infinitesimal external field to
break the symmetry and investigate the behavior in the limit of infinite
volume. We will discuss this issue in detail in Sec. III. Moreover, in the
case of $U(1)$ lattice gauge theory, the particle number probability
distribution function is zero except when the net particle number is zero
because of the $U(1)$ center symmetry. This is reasonable in the confinement
phase (disordered phase), but in the deconfinement phase (ordered phase), the
net particle number can be nonzero.

This paper is organized as follows. In the next section, we introduce the
probability distribution function of the particle density, which is given by
the canonical partition function. Then, we discuss problems caused by the
center symmetry in the calculation of the canonical partition function in Sec.
III. A method to calculate the canonical partition function by a saddle point
approximation is explained in Sec. IV. For the case of $U(1)$ lattice gauge
theory, the sign problem in the canonical approach is avoided at the same
time. We demonstrate this method for the $U(1)$ gauge theory with heavy
dynamical fermions. The results of the probability distribution function of
the particle density are shown in Sec. V. Then, the application to QCD is
discussed in Sec. VI. Our conclusions are given in Sec. VII.

##  II Canonical partition function and particle density probability
distribution function

The relation between the grand partition function $Z_{GC}(T,\mu)$ with the
chemical potential $\mu$ and the canonical partition function $Z_{C}(T,N)$
with the particle number $N$ is given by the fugacity expansion,

$\displaystyle Z_{GC}(T,\mu)=\sum_{N=-\infty}^{\infty}Z_{C}(T,N)e^{N\mu/T}$
(1)

at each temperature $T$. The left-hand side of this equation, $Z_{GC}(T,\mu)$,
is the normalization factor of the Boltzmann weight and is classified by $N$
in the right-hand side. Hence, $Z_{C}(T,N)e^{N\mu/T}$ can be regarded as a
weight factor for each $N$, and the probability distribution of $N$ is in
proportion to $Z_{C}(T,N)e^{N\mu/T}$.

Here, we introduce an effective potential as a function of the particle number
$N$ by the canonical partition function $Z_{C}(T,N)$,

$\displaystyle Z_{GC}(T,\mu)$ $\displaystyle=$
$\displaystyle\sum_{N=-\infty}^{\infty}e^{-V_{\rm eff}(T,N)},$ (2)
$\displaystyle V_{\rm eff}(T,N)$ $\displaystyle=$ $\displaystyle-\ln
Z_{C}(T,N)-\frac{N\mu}{T}.$ (3)

We note that there is an uncertainty in adding a constant to $V_{\rm eff}$,
which corresponds to an indefiniteness of a constant multiple of $Z_{GC}$. The
generation probability is maximized when the number of particles $N$ where
$V_{\rm eff}(T,N)$ is the minimum. We introduce the particle density $\rho$ as
$\rho=N/V$, where $V$ is the spatial volume. For sufficient large $V$, the
minimum point is obtained from the derivative of $V_{\rm eff}(T,V\rho)$,

$\displaystyle\frac{1}{V}\frac{\partial V_{\rm
eff}(T,V\rho)}{\partial\rho}=-\frac{1}{V}\frac{\partial\ln
Z_{C}(T,V\rho)}{\partial\rho}-\frac{\mu}{T}=0.$ (4)

Hence, the derivative $-(1/V)\partial\ln Z_{C}/\partial\rho$ is the $\mu/T$
where the particle density with the maximum generation probability is $\rho$.

The grand partition function of $N_{\rm f}$ flavor QCD is defined by

$\displaystyle Z_{GC}(T,\mu_{1},\mu_{2},\cdots,\mu_{N_{\rm f}})=\int{\cal
D}U\prod_{f=1}^{N_{\rm f}}\det M(\kappa_{f},\mu_{f}/T)\ e^{-S_{g}},$ (5)

where $\kappa_{f}$ is the hopping parameter and $\mu_{f}$ is the chemical
potential for the $f$th flavor. The action we study consists of the gauge
action and the fermion action,

$\displaystyle S=S_{g}+S_{q},$ (6)

where the gauge action is the standard plaquette action given by

$\displaystyle S_{g}=-6N_{\rm site}\,\beta\,\hat{P}$ (7)

with $N_{\rm site}=N_{s}^{3}\times N_{t}$ the space-time lattice volume and
$\hat{P}$ the plaquette operator. For the case of the standard plaquette gauge
action of $SU(N_{c})$ gauge theory, the gauge coupling parameter is
$\beta=2N_{c}/g^{2}$ and

$\displaystyle\hat{P}=\frac{1}{6N_{c}N_{\rm site}}\sum_{x,\,\mu<\nu}{\rm Re\
tr}\left[U_{x,\mu}U_{x+\hat{\mu},\nu}U^{\dagger}_{x+\hat{\nu},\mu}U^{\dagger}_{x,\nu}\right].$
(8)

Here, $U_{x,\nu}$ is the link variable in the $\nu$ direction at site $x$ and
$x+\hat{\nu}$ the next site in the $\nu$ direction from $x$. For fermions, we
adopt the standard Wilson fermion action given by

$\displaystyle S_{q}=\sum_{f=1}^{N_{\rm
f}}\sum_{x,\,y}\bar{\psi}_{x}^{(f)}\,M_{xy}(\kappa_{f},\mu_{f}/T)\,\psi_{y}^{(f)},$
(9)

where $\psi^{(f)}$ is the fermion field and $M_{xy}$ is the Wilson fermion
kernel

$\displaystyle M_{xy}(\kappa_{f},\mu_{f}/T)$ $\displaystyle=$
$\displaystyle\delta_{xy}-\kappa_{f}\left\\{\sum_{i=1}^{3}\left[(1-\gamma_{i})\,U_{x,i}\,\delta_{y,x+\hat{i}}+(1+\gamma_{i})\,U_{y,i}^{\dagger}\,\delta_{y,x-\hat{i}}\right]\right.$
(10)
$\displaystyle\left.+e^{\mu_{f}a}(1-\gamma_{4})\,U_{x,4}\,\delta_{y,x+\hat{4}}+e^{-\mu_{f}a}(1+\gamma_{4})\,U_{y,4}^{\dagger}\,\delta_{y,x-\hat{4}}\right\\}.$

The lattice spacing is $a$. Since $T=(N_{t}a)^{-1}$, the dimensionless ratio
$\mu_{f}/T=\mu_{f}aN_{t}$.

To discuss the chemical potential dependence of the fermion determinant, we
perform hopping parameter expansion in the vicinity of the heavy fermion
limit, $\kappa=0$. (See e.g. Sec. 5.1.3 of Ref. Montvay94  or Sec. 11 of Ref.
Rothe92 .) For each flavor, we have

$\displaystyle\ln\left[\frac{\det M(\kappa,\mu/T)}{\det
M(0,0)}\right]\;=\;\sum_{n=1}^{\infty}\frac{1}{n!}\left[\frac{\partial^{n}\ln\det
M}{\partial\kappa^{n}}\right]_{\kappa=0}\kappa^{n}\;=\;\sum_{n=1}^{\infty}\frac{B_{n}}{n!}\,\kappa^{n},$
(11)

with $\det M(0,0)=1$ and

$\displaystyle B_{n}$ $\displaystyle\equiv$
$\displaystyle\left[\frac{\partial^{n}\ln\det
M}{\partial\kappa^{n}}\right]_{\kappa=0}=(-1)^{n+1}(n-1)!\,{\rm
tr}\left[\left(M^{-1}\,\frac{\partial
M}{\partial\kappa}\right)^{n}\right]_{\kappa=0}$ (12) $\displaystyle=$
$\displaystyle(-1)^{n+1}(n-1)!\,{\rm tr}\left[\left(\frac{\partial
M}{\partial\kappa}\right)^{n}\right],$

where $(\partial M/\partial\kappa)_{xy}$ is the hopping term following
$\kappa_{f}$ in the right-hand side of Eq. (10). Nonvanishing contributions to
the trace of Eq. (12) appear only when the product of the hopping terms forms
a connected closed loop in the space-time. Therefore, $B_{n}$ are the sum of
connected $n$-step Wilson loops Saito:2011fs . We classify them into Wilson-
loop-type terms and Polyakov-loop-type terms. The former are independent of
the boundary conditions, and the latter are closed by the antiperiodic
boundary condition for the time direction. The leading-order (LO) contribution
consists of the smallest Wilson-loop-type term, plaquette $\hat{P}$, defined
by Eq. (8), and the smallest Polyakov-loop-type term, Polyakov loop
$\hat{\Omega}$, defined by

$\displaystyle\hat{\Omega}=\frac{1}{N_{c}N_{s}^{3}}\sum_{\vec{x}}{\rm
tr}\left[U_{\vec{x},4}U_{\vec{x}+\hat{4},4}U_{\vec{x}+2\cdot\hat{4},4}\cdots
U_{\vec{x}+(N_{t}-1)\cdot\hat{4},4}\right],$ (13)

where $\sum_{\vec{x}}$ means a summation over one time slice.

Because of the antiperiodic boundary condition and gamma matrices in the
hopping terms, up to the LO contributions for $SU(N_{c})$ gauge theory
Saito:2013vja , Eq. (11) reads

$\displaystyle\ln\left[\prod_{f=1}^{N_{\rm f}}\det
M(\kappa_{f},\mu_{f}/T)\right]$ $\displaystyle=$ $\displaystyle 96N_{c}N_{\rm
site}\sum_{f=1}^{N_{\rm f}}\kappa_{f}^{4}\hat{P}$ (14)
$\displaystyle+2^{N_{t}+1}N_{c}N_{s}^{3}\left(\sum_{f=1}^{N_{\rm
f}}\kappa_{f}^{N_{t}}e^{\mu_{f}/T}\hat{\Omega}+\sum_{f=1}^{N_{\rm
f}}\kappa_{f}^{N_{t}}e^{\mu_{f}/T}\hat{\Omega}^{*}\right)+\cdots.$

The first term that is proportional to $\hat{P}$ can be absorbed into the
gauge action by a shift $\beta\rightarrow\beta^{*}$ with

$\displaystyle\beta^{*}=\beta+16N_{c}\sum_{f=1}^{N_{\rm f}}\kappa_{f}^{4}.$
(15)

The next-leading order (NLO) contributions in the Wilson-loop-type terms are
the six-step Wilson loops and are typical operators in improved gauge actions.
Thus, the contributions of these NLO operators can also be absorbed by a shift
of improvement parameters of the gauge action. Because a shift in improvement
parameters only affects the lattice discretization errors, the six-step Wilson
loop terms will not affect characteristic features of the system in the
continuum limit. Moreover, the coefficients of the NLO terms are very small
compared to improvement parameters of typical improved actions
Wakabayashi:2021eye . In contrast, the terms proportional to the Polyakov-
loop-type terms act like external magnetic fields in spin models and thus may
change the nature of the phase transition. The higher-order contribution is
discussed in Refs. Ejiri:2019csa ; Kiyohara:2021smr ; Wakabayashi:2021eye . We
moreover assume that $N_{t}$ is an even number in this study.

As seen in Eq. (14), the Wilson-loop-type terms do not depend on the chemical
potential. Since each term of $U_{x,4}$ in Eq. (10) always appears as a
combination of $e^{\mu_{f}a}U_{x,4}$ or $e^{-\mu_{f}a}U^{{\dagger}}_{x,4}$,
the chemical potential dependence of the expansion term containing $n_{+}$
$U_{x,4}$ and $n_{-}$ $U^{{\dagger}}_{x,4}$ is $e^{\mu_{f}a(n_{+}-n_{-})}$.
Because $n_{+}$ and $n_{-}$ satisfy $n_{+}-n_{-}=N_{t}N$ for a closed loop
with a winding number $N$, Wilson-loop-type terms do not depend on the
chemical potential, and Polyakov-loop-type terms closed by the antiperiodic
boundary condition with the winding number $N$ are in proportion to
$e^{NN_{t}\mu_{f}a}=e^{N\mu_{f}/T}$. Thus, the expansion terms can be
classified by the winding number $N$,

$\displaystyle\ln\left[\det
M(\kappa,\mu_{f}/T)\right]=\sum_{N=-\infty}^{\infty}C_{N}\,e^{N\mu_{f}/T}.$
(16)

$C_{N}e^{N\mu_{f}/T}$ is the sum of the terms whose winding number is $N$ in
the hopping parameter expansion. If we classify the expansion terms by the
winding number, the canonical partition function is obtained. For example,
when $N_{\rm f}=2$ with the same $\mu$ and $\kappa$, $Z_{C}(T,N)$ can be
calculated from the equation,

$\displaystyle Z_{GC}(T,\mu)$ $\displaystyle=$
$\displaystyle\int\mathcal{D}U(\mathrm{det}M(\kappa,\mu/T))^{2}e^{-S_{g}}=\int{\cal
D}U\exp\left[\sum_{m=-\infty}^{\infty}C_{m}e^{m\mu/T}+\sum_{n=-\infty}^{\infty}C_{n}e^{n\mu/T}\right]e^{-S_{g}}$
(17) $\displaystyle=\int{\cal
D}U\sum_{l=0}^{\infty}\frac{1}{l!}\left[\sum_{m=-\infty}^{\infty}C_{m}e^{m\mu/T}+\sum_{n=-\infty}^{\infty}C_{n}e^{n\mu/T}\right]^{l}e^{-S_{g}}$
$\displaystyle\equiv\sum_{N=-\infty}^{\infty}\left[\int{\cal
D}U\,{\det}_{N}M\,e^{-S_{g}}\right]e^{N\mu/T}=\sum_{N=-\infty}^{\infty}Z_{C}(T,N)\,e^{N\mu/T}$

for each term of $e^{N\mu/T}$, if we perform path integral over gauge
configurations. Here, $\det_{N}M$ means the expansion term of $[\det
M(\kappa,\mu/T)]^{N_{\rm f}}$ with respect to $N$. Therefore, the fugacity
expansion is the winding number $N$ expansion. Then, the chemical potential
enlarges the contributions having positive winding number and suppress the
contributions having negative winding number.

##  III Center symmetry and the canonical partition function

###  III.1 $SU(3)$ gauge theory

The quenched QCD, in which no dynamical fermions are included, has $Z_{3}$
center symmetry. The centers of $SU(3)$ group are $\\{I,\omega I$,
$\omega^{2}I\\}$, where $\omega=e^{2\pi i/3}$ and $I$ is the $3\times 3$ unit
matrix. The center transformation is defined as

$\displaystyle U_{(\vec{x},t),4}\to\omega\ U_{(\vec{x},t),4}$ (18)

for all $\vec{x}$ in one time slice $t$. Under the center transformation, the
system is symmetric, i.e. the gauge action and integral measure ${\cal D}U$
are invariant. However, the expectation value of Polyakov loop $\hat{\Omega}$
changes as

$\displaystyle\langle\hat{\Omega}\rangle\to\omega\langle\hat{\Omega}\rangle.$
(19)

If this symmetry is maintained, any expectation values do not change under the
center transformation. Then,
$\langle\hat{\Omega}\rangle=\omega\langle\hat{\Omega}\rangle$ must be
satisfied. Thus, $\langle\hat{\Omega}\rangle$ is zero. Therefore, the Polyakov
loop is regarded as an order parameter of the spontaneous breaking of the
center symmetry.

Similarly, under the $Z_{3}$ center transformation, the canonical partition
function of the particle number $N$, $Z_{C}(T,N)$, changes as

$\displaystyle Z_{C}(T,N)\to\omega^{N}Z_{C}(T,N)=e^{2\pi Ni/3}Z_{C}(T,N),$
(20)

since $Z_{C}(T,N)$ is given as the sum of the expectation values of Polyakov-
loop-type operator with a winding number of $N$ as discussed in the previous
section. Because the grand partition function is given by the fugacity
expansion, $Z_{GC}(T,\mu)$ changes as

$\displaystyle\sum_{N=-\infty}^{\infty}Z_{C}(T,N)e^{N\mu/T}\to\sum_{N=-\infty}^{\infty}Z_{C}(T,N)e^{N\mu/T+2\pi
Ni/3},$ (21)

and is not invariant under the center transformation when including dynamical
fermions. At the same time, this implies

$\displaystyle Z_{GC}(T,\mu)=Z_{GC}(T,\mu+2\pi iT/3),$ (22)

since $S_{g}$ and the integral measure do not change under the center
transformation, where we extend the real $\mu$ to a complex number. This
symmetry is called Roberge-Weiss symmetry Roberge:1986mm .

Moreover, this indicates $Z_{C}(T,N)=0$ expect when $N$ is a multiple of 3.
Since $\omega^{3}=1$ and $1+\omega+\omega^{2}=0$,

$\displaystyle Z_{GC}(T,\mu)$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left(Z_{GC}(T,\mu)+Z_{GC}(T,\mu+2\pi
iT/3)+Z_{GC}(T,\mu+4\pi iT/3)\right)$ (23) $\displaystyle=$
$\displaystyle\sum_{n=-\infty}^{\infty}Z_{C}(T,3n)e^{3n\mu/T},$

where $n$ is an integer. This means that the number of particles can only
exist in multiples of three. This is reasonable in the confinement phase where
only baryons and mesons exist. However, in the deconfinement phase, there
should be states of all particle numbers. The order parameter of spontaneous
symmetry breaking has a nonzero expectation value when one vacuum is selected
from multiple vacua in the limit of infinite volume. Thus, the center symmetry
is never broken in an actual simulation with finite volume.

###  III.2 $U(1)$ gauge theory

This problem is more serious in $U(1)$ gauge theory. For the case of $U(1)$
gauge theory, the centers are all elements of $U(1)$ group because of an
Abelian group. The link variable is given as $U_{x,\nu}=e^{i\theta_{x,\nu}}$,
where $\theta_{x,\nu}$ is a real variable that takes the value of
$-\pi<\theta_{x,\nu}\leq\pi$ on a link $(x,\nu)$. The gauge action is given by
$S_{g}=-6N_{\rm site}\,\beta\,\hat{P},$ with $\beta=1/g^{2}$ the gauge
coupling parameter and $\hat{P}$ the plaquette operator,

$\displaystyle\hat{P}=\frac{1}{6N_{\rm
site}}\sum_{x,\,\mu<\nu}\cos\left[\theta_{x,\mu}+\theta_{x+\hat{\mu},\nu}-\theta_{x+\hat{\nu},\mu}-\theta_{x,\nu}\right]$
(24)

for the case of the standard plaquette action. When we perform the hopping
parameter expansion, the LO contribution of the fermion determinant is given
by

$\displaystyle\ln\left[\prod_{f=1}^{N_{\rm f}}\det
M(\kappa_{f},\mu_{f}/T)\right]$ $\displaystyle=$ $\displaystyle 96N_{\rm
site}\sum_{f=1}^{N_{\rm f}}\kappa_{f}^{4}\hat{P}$ (25)
$\displaystyle+2^{N_{t}+1}N_{s}^{3}\left(\sum_{f=1}^{N_{\rm
f}}\kappa_{f}^{N_{t}}e^{\mu_{f}/T}\hat{\Omega}+i\sum_{f=1}^{N_{\rm
f}}\kappa_{f}^{N_{t}}e^{\mu_{f}/T}\hat{\Omega}^{*}\right)+\cdots.$

The Polyakov loop for $U(1)$ gauge theory is

$\displaystyle\hat{\Omega}=\frac{1}{N_{s}^{3}}\sum_{\vec{x}}\exp\left[i(\theta_{\vec{x},4}+\theta_{\vec{x}+\hat{4},4}+\theta_{\vec{x}+2\hat{4},4}+\cdots+\theta_{\vec{x}+N_{t}\hat{4},4})\right].$
(26)

Under the $U(1)$ center transformation: $U_{(\vec{x},t),4}\to
e^{i\eta}U_{(\vec{x},t),4}$ on one time slice $t$, $Z_{C}(T,N)$ changes as

$\displaystyle Z_{C}(T,N)\to e^{iN\eta}Z_{C}(T,N).$ (27)

Since $S_{g}$ and the integral measure are invariant,
$Z_{C}(T,N)=e^{iN\eta}Z_{C}(T,N)$. Thus, the canonical partition function is

$\displaystyle
Z_{C}(T,N)=\frac{1}{2\pi}\int_{0}^{2\pi}e^{iN\eta}Z_{C}(T,N)\,d\eta=0,$ (28)

except for $N=0$. This means that the existence probability of particles that
interact with the gauge field is zero. Moreover, the grand partition function
is

$\displaystyle Z_{GC}(T,\mu_{q})=Z_{C}(T,0)$ (29)

and does not depend on the chemical potential due to the $U(1)$ center
symmetry. Since the center symmetry cannot be broken in numerical simulations
performed with finite volume, the canonical partition function cannot be
calculated correctly.

###  III.3 Polyakov loop in $U(1)$ gauge theory

Figure 1: Distribution of Polyakov loops in the complex plane generated by
simulations of $U(1)$ lattice gauge theory at $\beta=0.90$ (left) and 1.10
(right). The temporal lattice size $N_{t}$ is 4. Purple, blue, green and red
symbols are the results of the spatial lattice size $N_{s}=12$, 16, 24, and
32, respectively.

This problem is essentially the same as the problem that the expectation value
of the Polyakov loop is always zero in a practical simulation with finite
volume. In the case of heavy fermions, i.e. $\kappa$ is sufficiently small,
the canonical partition function with $N=1$ is in proportion to
$\langle\hat{\Omega}\rangle$, since the leading contribution is the Polyakov
loop term in the hopping parameter expansion of $\ln\det M$ with the winding
number one, and $\det_{1}M$ in the equation of $Z_{C}(T,1)$, e.g. Eq. (17), is
proportional to $\hat{\Omega}$. In Fig. 1, we plot complex values of the
Polyakov loop on each configuration. The configurations are generated by a
heat bath method of $U(1)$ lattice gauge theory Creutz:1983ev  with the
standard plaquette gauge action and no dynamical fermions. The left panel
shows the result of the confinement (symmetric) phase at $\beta=0.90$, and the
right panel shows the result of the deconfinement (broken) phase at
$\beta=1.10$. The temporal lattice size $N_{t}$ is fixed to be 4. We adopt
four spatial lattice sizes $N_{s}$. Purple, blue, green and red symbols are
the results of $N_{s}=12$, 16, 24 and 32, respectively. In the confinement
phase, the Polyakov loop distribute around $\Omega=0$. As the volume
increases, the width of the distribution decreases. In the deconfinement
phase, the Polyakov loop is distributed on a circle and there is no volume
dependence in the distribution.

The Polyakov loop is an order parameter of the deconfinement phase transition.
In the deconfinement phase, the center symmetry is spontaneously broken and
$\langle\hat{\Omega}\rangle$ should be nonzero. However, because of the center
symmetry, the probability distribution is symmetric under the $U(1)$
transformation, $\hat{\Omega}\to e^{i\eta}\hat{\Omega}$ for an arbitrary real
number $\eta$, as shown in Fig. 1. Therefore, the expectation value of
$\hat{\Omega}$ is always zero even in the broken phase. Because the symmetry
is not broken in an simulation, to discuss spontaneous symmetry breaking, it
is required to break the center symmetry adding an explicit breaking term in
the action. Then, the breaking-term dependence and spatial-volume dependence
are investigated. If $\langle\hat{\Omega}\rangle$ is nonzero in the double
limit of zero breaking term and infinite volume, we identify that spontaneous
symmetry breaking has occurred.

To break the center symmetry, we add one flavor of heavy dynamical fermion. We
approximate the fermion determinant with the leading-order term of the hopping
parameter $\kappa$ expansion given in Eq. (14), since we investigate the limit
of $\kappa\to 0$. The plaquette term can be absorbed into the gauge action by
shifting $\beta\to\beta^{*}=\beta+16\kappa^{4}$. The expectation value is
computed by the reweighting method with the reweighting factor including the
Polyakov-loop term,

$\displaystyle\langle{\rm Re}\hat{\Omega}\rangle_{(\beta,\kappa)}$
$\displaystyle=$ $\displaystyle\frac{1}{Z_{GC}}\int{\cal D}U\,{\rm
Re}\hat{\Omega}\,\det M\,e^{-S_{g}(\beta)}$ (30) $\displaystyle\approx$
$\displaystyle\frac{1}{Z_{GC}}\int{\cal D}U\,{\rm Re}\hat{\Omega}\,e^{\epsilon
V{\rm Re}\hat{\Omega}}\,e^{-S_{g}(\beta^{*})}=\frac{\left\langle{\rm
Re}\hat{\Omega}\,e^{\epsilon V{\rm
Re}\hat{\Omega}}\right\rangle_{(\beta^{*},0)}}{\left\langle e^{\epsilon V{\rm
Re}\hat{\Omega}}\right\rangle_{(\beta^{*},0)}},$

where $\epsilon=4\times 2^{N_{t}}\kappa^{N_{t}}$, $V=N_{s}^{3}$ and
$\langle\cdots\rangle_{(\beta^{*},0)}$ means the average over quenched
configurations at $\beta^{*}$. Using the data of $N_{t}=4$ in Fig. 1, the
expectation value of ${\rm Re}\hat{\Omega}$ is computed. The number of updates
is 1,000,000 for each $\beta$. The jackknife error is evaluated adopting an
appropriate bin size. We plot the results in Fig. 2 as a function of
$\kappa^{N_{t}}$. The left figure is the result of the confinement (symmetric)
phase at $\beta=0.90$. The results of $N_{s}=12$ (purple), 16 (blue), 24
(green) and 32 (red) are plotted. No volume dependence is observed for
$\langle{\rm Re}\hat{\Omega}\rangle$. Therefore, it is not necessary to take
the volume infinity limit, and $\langle{\rm
Re}\hat{\Omega}\rangle\sim\kappa^{N_{t}}$ in the confinement phase (symmetric
phase). Then, $\langle{\rm Re}\hat{\Omega}\rangle=0$ in the limit of
$\kappa\to 0$. On the other hand, the right figure is the result of the
deconfinement (broken) phase at $\beta=1.10$. In the deconfinement phase, the
Polyakov loop behaves as $\langle{\rm Re}\hat{\Omega}\rangle\sim
V\kappa^{N_{t}}$. Of course, if $V$ is fixed and the limit of $\kappa\to 0$ is
taken, $\langle{\rm Re}\hat{\Omega}\rangle$ becomes always zero. However, in
the thermodynamic limit, i.e. $V\to\infty$, and $\kappa\to 0$, $\langle{\rm
Re}\hat{\Omega}\rangle$ is expected to become a nontrivial finite value.111 In
this analysis, we ignored the difference between $\beta$ and $\beta^{*}$
because the difference is very small. To be precise, Figure 2 is the results
of $\beta^{*}=0.90$ and 1.10. Although $\langle{\rm Re}\hat{\Omega}\rangle$ at
$\kappa=0$ is zero within the statistical error due to the center symmetry, in
order to make it exactly zero at $\kappa=0$, we assumed that the generation
probability of $\hat{\Omega}$ and $-\hat{\Omega}$ is the same in the
calculation of the expectation value.

Figure 2: Expectation value of the Polyakov loop as a function of the hopping
parameter $\kappa$ measured at $\beta=0.90$ (left) and 1.10 (right) on
$N_{s}^{3}\times 4$ lattices with $N_{s}=12$ (purple), 16 (blue), 24 (green),
and 32 (red).

Furthermore, when $\kappa$ is sufficiently small, the expectation value of the
Polyakov loop can be evaluated by a Taylor expansion of $\kappa^{N_{t}}$
assuming the $U(1)$ center symmetry. Using the distribution function of the
Polyakov loop in the complex plane, the expectation value can be calculated as
follows,

$\displaystyle\langle{\rm Re}\hat{\Omega}\rangle_{(\beta,\kappa)}$
$\displaystyle=$ $\displaystyle\frac{1}{Z_{GC}}\int{\cal D}U\,{\rm
Re}\hat{\Omega}\,e^{\epsilon V{\rm
Re}\hat{\Omega}}\,e^{-S_{g}}=\frac{1}{2\pi}\iint|\Omega|\cos\phi\,e^{\epsilon
V|\Omega|\cos\phi}\,W(|\Omega|)\,d\phi\,d|\Omega|$ (31) $\displaystyle=$
$\displaystyle\frac{\epsilon
V}{2}\int|\Omega|^{2}\,W(|\Omega|)\,d|\Omega|+\cdots.$

Here, $W(|\Omega|)$ is the probability distribution when the Polyakov loop is
$|\Omega|e^{i\phi}$ in the simulation without dynamical fermions. This
probability distribution function is $U(1)$ symmetric, i.e. independent of
$\phi$, and is a function of the absolute value $|\Omega|$, as seen in Fig. 1.
The distribution function is normalized as
$\int_{0}^{\infty}W(|\Omega|)\,d|\Omega|=1$. The value of $\langle{\rm
Re}\hat{\Omega}^{n}\rangle$ can be also computed,

$\displaystyle\langle{\rm Re}\hat{\Omega}^{n}\rangle_{(\beta,\kappa)}$
$\displaystyle=$ $\displaystyle\frac{1}{Z_{GC}}\int{\cal D}U\,{\rm
Re}\hat{\Omega}^{n}\,e^{\epsilon V{\rm
Re}\hat{\Omega}}e^{-S_{g}}=\frac{1}{2\pi}\iint|\Omega|^{n}\cos(n\phi)\,e^{\epsilon
V|\Omega|\cos\phi}W(|\Omega|)\,d\phi\,d|\Omega|$ (32) $\displaystyle=$
$\displaystyle\frac{(\epsilon V)^{n}}{n!\
2^{n}}\int|\Omega|^{2n}\,W(|\Omega|)\,d|\Omega|+\cdots.$

The complex phase can be removed in the above equations. In the double limit
of $V\to\infty$, and $\kappa\to 0$, the value of $\epsilon V$ cannot be
determined. However, for example, the leading term of the ratio $\langle{\rm
Re}\hat{\Omega}^{4}\rangle/\langle{\rm Re}\hat{\Omega}^{2}\rangle^{2}$ does
not depend on $\epsilon V$. The volume dependence of such a quantity is
expected to small. The explicit breaking term is also important for the
calculation of the canonical partition function. As discussed below, this
method can be applied in computing the derivative of canonical partition
function.

##  IV Canonical partition function with a saddle point approximation

We calculate the canonical partition function following Ref. Ejiri:2008xt .
Because the fugacity expansion is basically a Laplace transform, the canonical
partition function is obtained by the inverse Laplace transform,

$\displaystyle{\cal Z}_{\rm
C}(T,N)=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{-N(\mu_{0}/T+i\mu_{I}/T)}{\cal
Z}_{\rm GC}(T,\mu_{0}+i\mu_{I})\ d\left(\frac{\mu_{I}}{T}\right),$ (33)

where $\mu_{0}$ is arbitrary complex number. $\mu_{I}/T$ is an integral
variable introduced to constrain the number of particles $N$. Regarding as the
imaginary part of $\mu/T$, we extend the $\mu/T$ in Eq. (1) to a complex
value. This equation is equivalent to performing the complex integral on the
complex plane of $z$ along a path parallel to the imaginary axis,

$\displaystyle Z_{C}(T,N)=\frac{1}{2\pi}\int Z_{GC}(T,zT)e^{-Nz}dz,$ (34)

where $z=(\mu_{0}+i\mu_{I})/T$.

We consider the case of degenerate $N_{\rm f}$ flavor QCD for simplicity. The
grand partition function:

$\displaystyle Z_{GC}(T,zT)$
$\displaystyle=\int\mathcal{D}U(\mathrm{det}M(\kappa,z))^{N_{\rm
f}}e^{-S_{g}}$ (35)

can not be computed by Monte Carlo simulations. Thus, we calculate the
following ratio,

$\displaystyle\frac{{\cal Z}_{\rm GC}(T,\mu)}{{\cal Z}_{\rm GC}(T,0)}$
$\displaystyle=$ $\displaystyle\frac{1}{{\cal Z}_{\rm GC}}\int{\cal
D}U\left(\frac{\det M(\kappa,\mu/T)}{\det M(\kappa_{0},0)}\right)^{N_{\rm
f}}(\det M(\kappa_{0},0))^{N_{\rm f}}e^{-S_{g}}$ (36) $\displaystyle=$
$\displaystyle\left\langle\left(\frac{\det M(\kappa,\mu/T)}{\det
M(\kappa_{0},0)}\right)^{N_{\rm f}}\right\rangle_{(T,\mu=0,\kappa_{0})}.$

$Z_{\rm GC}(T,0)$ is the grand partition function at the simulation point of
$\mu=0$. $\langle\cdots\rangle_{(T,\mu=0,\kappa_{0})}$ is the expectation
value at $\mu=0$ with the hopping parameter $\kappa_{0}$. Substituting this
equation into Eq. (34), the canonical partition function becomes

$\displaystyle Z_{C}(T,N)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}Z_{GC}(T,0)\int\left\langle
e^{VD(z)}e^{-Nz}\right\rangle_{(T,\mu=0)}dz$ (37) $\displaystyle=$
$\displaystyle\frac{1}{2\pi}Z_{GC}(T,0)\left\langle\int
e^{V(D(z)-\frac{N}{V}z)}dz\right\rangle_{(T,\mu=0)},$

where $V=N_{s}^{3}$ and

$\displaystyle e^{VD(z)}\equiv\left(\frac{\det M(\kappa,z)}{\det
M(\kappa_{0},0)}\right)^{N_{\rm f}}.$ (38)

In the second line of Eq. (37), we changed the calculation order of the
complex integral of $z$ and the path integral of the gauge field.

Here, the particle number density is expressed as $\rho=N/V$ in lattice units.
In physical units, the density $\rho_{\rm phy}$ is given by $\rho_{\rm
phy}/T^{3}=\rho N_{t}^{3}$. For simplicity, we integrate $z$ by a saddle point
approximation Ejiri:2008xt . This approximation is valid for large volumes.
The saddle point condition is

$\displaystyle\frac{d}{dz}\left[D(z)-\rho z\right]=\frac{dD(z)}{dz}-\rho=0.$
(39)

Let $z_{0}=x_{0}+iy_{0}$ be the saddle point where $z$ satisfies the saddle
point condition.

$\displaystyle{\cal Z}_{\rm C}(T,\rho V)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}{\cal Z}_{\rm
GC}(T,0)\left\langle\int_{-\pi}^{\pi}e^{-N(z_{0}+ix)}\left(\frac{\det
M(\kappa,z_{0}+ix)}{\det M(\kappa_{0},0)}\right)^{N_{\rm
f}}dx\right\rangle_{(T,\mu=0,\kappa_{0})}$ (40) $\displaystyle=$
$\displaystyle\frac{1}{2\pi}{\cal Z}_{\rm
GC}(T,0)\left\langle\int_{-\pi}^{\pi}\exp\left[V\left(D(z_{0})-\rho
z_{0}-\frac{1}{2}D^{\prime\prime}(z_{0})x^{2}+\cdots\right)\right]dx\right\rangle_{(T,\mu=0,\kappa_{0})}$
$\displaystyle\approx$ $\displaystyle\frac{1}{\sqrt{2\pi}}{\cal Z}_{\rm
GC}(T,0)\left\langle\exp\left[V\left(D(z_{0})-\rho
z_{0}\right)\right]e^{-i\alpha/2}\sqrt{\frac{1}{V|D^{\prime\prime}(z_{0})|}}\right\rangle_{(T,\mu=0,\kappa_{0})},$

where $D^{\prime\prime}(z)=d^{2}D(z)/dz^{2}$ and
$D^{\prime\prime}(z)=|D^{\prime\prime}(z)|e^{i\alpha}$. The derivative with
respect to $\rho$ is obtained by the following equation,

$\displaystyle-\frac{1}{V}\frac{\partial\ln{\cal Z}_{\rm C}(T,\rho
V)}{\partial\rho}\approx\frac{\left\langle z_{0}\
\exp\left[V\left(D(z_{0})-\rho
z_{0}\right)\right]e^{-i\alpha/2}\sqrt{\frac{1}{V|D^{\prime\prime}(z_{0})|}}\right\rangle_{(T,\mu=0,\kappa_{0})}}{\left\langle\exp\left[V\left(D(z_{0})-\rho
z_{0}\right)\right]e^{-i\alpha/2}\sqrt{\frac{1}{V|D^{\prime\prime}(z_{0})|}}\right\rangle_{(T,\mu=0,\kappa_{0})}}.$
(41)

This equation can be calculated by the Monte Carlo method.

##  V $U(1)$ lattice gauge theory with heavy fermions

###  V.1 Canonical partition function of heavy fermions

We demonstrate the calculation of the canonical partition function in $U(1)$
lattice gauge theory with degenerate $N_{\rm f}$ flavors of heavy dynamical
fermions. $D(z)$ in Eq. (37) is evaluated by the hopping parameter expansion,

$\displaystyle D(z)=96N_{t}N_{\rm f}\kappa^{4}\hat{P}+2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left[e^{z}\hat{\Omega}+e^{-z}\hat{\Omega}^{{\dagger}}\right]+\cdots.$
(42)

Assuming that $\kappa$ is sufficiently small, we approximate $D(z)$ by the
terms $\hat{P}$ and $\hat{\Omega}$. The first and second derivatives of $D(z)$
with respect to $z$ are

$\displaystyle\frac{\partial D(z)}{\partial z}$ $\displaystyle=$
$\displaystyle 2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left[e^{z}\hat{\Omega}-e^{-z}\hat{\Omega}^{{\dagger}}\right]+\cdots,$
(43) $\displaystyle\frac{\partial^{2}D(z)}{\partial z^{2}}$ $\displaystyle=$
$\displaystyle 2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left[e^{z}\hat{\Omega}+e^{-z}\hat{\Omega}^{{\dagger}}\right]+\cdots.$
(44)

The saddle point $z_{0}=x_{0}+iy_{0}$ is defined as

$\displaystyle\left[\frac{\partial D(z)}{\partial z}-\rho\right]_{z=z_{0}}=0.$
(45)

Thus,

$\displaystyle 2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left(e^{z_{0}}\hat{\Omega}-e^{-z_{0}}\hat{\Omega}^{{\dagger}}\right)=\rho.$
(46)

Because $\rho$ is a real number, the complex phase of
$(e^{iy_{0}}\hat{\Omega}-e^{-iy_{0}}\hat{\Omega}^{{\dagger}})$ is zero. The
imaginary part of saddle point $y_{0}$ is determined as

$\displaystyle\tan
y_{0}=-\frac{\mathrm{Im}\hat{\Omega}}{\mathrm{Re}\hat{\Omega}}\hskip
14.22636pt{\rm or}\hskip
14.22636pty_{0}=-\arctan\left(\frac{\mathrm{Im}\hat{\Omega}}{\mathrm{Re}\hat{\Omega}}\right)=-\mathrm{Arg}\hat{\Omega}.$
(47)

Substituting $y_{0}$ into Eq. (46),

$\displaystyle 2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left[e^{x_{0}}|\hat{\Omega}|-e^{-x_{0}}|\hat{\Omega}|\right]=4\times
2^{N_{t}}N_{\rm f}\kappa^{N_{t}}|\hat{\Omega}|\sinh x_{0}=\rho.$ (48)

Solving $x_{0}$,

$\displaystyle x_{0}$ $\displaystyle=$
$\displaystyle\mathrm{arcsinh}\left(\frac{\rho}{4\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}|\hat{\Omega}|}\right)$ (49) $\displaystyle=$
$\displaystyle\ln\left(\frac{\rho}{4\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}|\hat{\Omega}|}+\sqrt{\left(\frac{\rho}{4\times
2^{N_{t}}N_{\rm f}\kappa^{N_{t}}|\hat{\Omega}|}\right)^{2}+1}\right).$

The second derivative of $D(z)$ at $z=z_{0}$, $D^{\prime\prime}(z_{0})$, is
given by

$\displaystyle D^{\prime\prime}(z_{0})=2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left[e^{x_{0}}|\hat{\Omega}|+e^{-x_{0}}|\hat{\Omega}|\right]=4\times
2^{N_{t}}N_{\rm f}\kappa^{N_{t}}|\hat{\Omega}|\cosh x_{0}.$ (50)

Similarly, $D(z)$ at $z=z_{0}$ is

$\displaystyle D(z_{0})=96N_{t}N_{\rm f}\kappa^{4}\hat{P}+4\times
2^{N_{t}}N_{\rm f}\kappa^{N_{t}}|\hat{\Omega}|\cosh x_{0}.$ (51)

In this approximation, $x_{0}$, $D(z_{0})$, and $D^{\prime\prime}(z_{0})$ are
real positive numbers. We substitute these equations into Eq. (41). The
derivative of $Z_{C}$ is

$\displaystyle-\frac{1}{V}\frac{\partial\mathrm{ln}Z_{C}(T,N)}{\partial\rho}\approx\frac{\left\langle
z_{0}\exp(F+i\theta)\right\rangle}{\left\langle\exp(F+i\theta)\right\rangle},$
(52)

where $F$ and $\theta$ are defined as

$\displaystyle F$ $\displaystyle=$ $\displaystyle V(D(z_{0})-\rho
x_{0})-\frac{1}{2}\ln[VD^{\prime\prime}(z_{0})],$ (53) $\displaystyle\theta$
$\displaystyle=$
$\displaystyle-V\rho\,\mathrm{Arg}\hat{\Omega}=-V\rho\,\arctan\left(\frac{{\rm
Im}\hat{\Omega}}{{\rm Re}\hat{\Omega}}\right).$ (54)

The Monte Carlo simulation is performed without dynamical fermions,
$\kappa_{0}=0$ and $\det M(0,0)=1$, and $\langle\cdots\rangle$ means the
expectation value of the quenched simulations. The term proportional to
$\hat{P}$ in Eq. (51) can be absorbed into the gauge action by shifting
$\beta\to\beta^{*}=\beta+16N_{\rm f}\kappa^{4}$ in the quenched simulation,
and the shift of $\beta$ is very small in this calculation. Therefore, we omit
Wilson-loop-type terms such as the plaquette term and concentrate on the
effect of the Polyakov-loop term. Then, $x_{0}$ and $F$ are real functions of
$\rho$ and the absolute value of $\Omega$. (See Eqs. (49), (50) and (51).) On
the other hand, $\theta$ is a real function of $\rho$ and the complex phase of
$\Omega$.

We calculate the derivative of the canonical partition function by classifying
the configurations by the value of $|\Omega|$ in the Monte Carlo simulation.

$\displaystyle-\frac{1}{V}\frac{\partial\ln Z_{C}(T,V\rho)}{\partial\rho}$
$\displaystyle\approx$ $\displaystyle\frac{\int\langle
z_{0}\,\exp[F+i\theta]\rangle_{|\Omega|}\,W(|\Omega|)\,d|\Omega|}{\int\langle\exp[F+i\theta]\rangle_{|\Omega|}\,W(|\Omega|)\,d|\Omega|}$
(55) $\displaystyle=$
$\displaystyle\frac{\int(x_{0}\langle\cos\theta\rangle_{|\Omega|}+\langle
y_{0}\sin\theta\rangle_{|\Omega|})\,e^{F}\,W(|\Omega|)\,d|\Omega|}{\int
e^{F}\langle\cos\theta\rangle_{|\Omega|}\,W(|\Omega|)\,d|\Omega|},$

where $\langle\cdots\rangle_{|\Omega|}$ means that each configuration is
classified by the value of $|\Omega|$ and the expectation value is calculated
for each value of $|\Omega|$, i.e.

$\displaystyle\langle\cdots\rangle_{|\Omega|=x}=\frac{\langle\cdots\delta(|\hat{\Omega}|-x)\rangle}{\langle\delta(|\hat{\Omega}|-x)\rangle},$
(56)

and $W(|\Omega|)$ is the probability distribution function of $|\Omega|$.
However, because $\theta=-Ny_{0}$ and
$y_{0}=-\mathrm{Arg}\hat{\Omega}\equiv-\phi$,
$\langle\cos\theta\rangle_{|\Omega|}=\langle\cos N\phi\rangle_{|\Omega|}=0$.
Namely, the phase average $\langle\cos\theta\rangle_{|\Omega|}$ is exactly
zero due to the $U(1)$ symmetric distribution function of $\Omega$. So to
speak, this is the ultimate sign problem. Therefore, $\partial\ln
Z_{C}/\partial\rho$ is indefinite.

Here, we neglect a term $\int e^{F}\langle
y_{0}\sin\theta\rangle_{|\Omega|}W(|\Omega|)\,d|\Omega|$ in Eq. (55). Because
of the $U(1)$ center symmetry, this can be computed as follows:

$\displaystyle\langle y_{0}\sin\theta\rangle_{|\Omega|}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{c}^{c+2\pi}\phi\sin(N\phi)\,d\phi=\frac{1}{2\pi}\left[-\frac{\phi}{N}\cos(N\phi)+\frac{1}{N^{2}}\sin(N\phi)\right]_{c}^{c+2\pi}$
(57) $\displaystyle=$ $\displaystyle-\frac{1}{N}\cos(Nc).$

This value of $\langle y_{0}\sin\theta\rangle_{|\Omega|}$ changes if the upper
and lower bounds of the integration are changed while the integration range
remains $2\pi$. We thus define this quantity as the product of an appropriate
convergence factor, which satisfies $\lim_{\phi\to\infty}e^{-\xi\phi}=0$.
Then, this term vanishes,

$\displaystyle\langle y_{0}\sin\theta\rangle_{|\Omega|}$ $\displaystyle\equiv$
$\displaystyle\lim_{\xi\to 0,\
m\to\infty}\frac{1}{2m\pi}\int_{0}^{2m\pi}\phi\sin(N\phi)\,e^{-\xi\phi}\,d\phi$
(58) $\displaystyle=$ $\displaystyle\lim_{\xi\to 0,\
m\to\infty}\frac{1}{2m\pi}\left[-\frac{e^{-\xi\phi}}{(N^{2}+\xi^{2})^{2}}\left\\{\left(N(N^{2}+\xi^{2})\phi+2N\xi\right)\cos(N\phi)\right.\right.$
$\displaystyle\hskip
56.9055pt\left.\left.+\left(\xi(N^{2}+\xi^{2})\phi-N^{2}+\xi^{2}\right)\sin(N\phi)\right\\}\right]_{0}^{2m\pi}$
$\displaystyle=$ $\displaystyle\lim_{\xi\to 0,\
m\to\infty}\left[\frac{-Ne^{-2m\pi\xi}}{N^{2}+\xi^{2}}+\frac{N\xi(1-e^{-2m\pi\xi})}{m\pi(N^{2}+\xi^{2})^{2}}\right]=0,$

where $\xi$ is a real number and $m$ is an integer number. Moreover, $\langle
y_{0}\sin\theta\cos^{n}\phi\rangle_{|\Omega|}$ also vanishes for any positive
integer $n$, since it is given by the sum of $\int\phi\sin(k\phi)\,d\phi$ with
an appropriate integer $k$. When we break the center symmetry adding a heavy
dynamical fermion, the calculation of terms like $\langle
y_{0}\sin\theta\cos^{n}\phi\rangle_{|\Omega|}$ is required, but such terms do
not contribute to Eq. (55). Therefore, we neglect the term of $\langle
y_{0}\sin\theta\rangle_{|\Omega|}$.

We use the method explained in Sec. III.3 to compute
$\langle\cos\theta\rangle_{|\Omega|}$ avoiding the sign problem. An additional
heavy fermion with a small hopping parameter $\kappa_{h}$ is introduced to
break the center symmetry. The fermion determinant is approximated by
$e^{\epsilon V{\rm Re}\hat{\Omega}}$, where $\hat{\Omega}$ is the Polyakov-
loop operator and $\epsilon=4\times 2^{N_{t}}\kappa_{h}^{N_{t}}$, and the
Wilson-loop-type terms are disregarded. When $\epsilon$ is small, the
denominator of Eq. (55) is computed as follows,

$\displaystyle\langle e^{F}\cos\theta\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{Z}\int{\cal D}Ue^{F}\cos\theta\,e^{\epsilon V{\rm
Re}\hat{\Omega}}e^{-S_{g}}=\frac{1}{Z}\int{\cal
D}Ue^{F}\cos(N\phi)\,e^{\epsilon V|\hat{\Omega}|\cos\phi}e^{-S_{g}}$ (59)
$\displaystyle=$ $\displaystyle\frac{1}{2\pi}\iint
e^{F}\cos(N\phi)\,e^{\epsilon V|\Omega|\cos\phi}W(|\Omega|)\,d\phi\,d|\Omega|$
$\displaystyle=$ $\displaystyle\int
e^{F}W(|\Omega|)\,\frac{1}{2^{N}N!}\,(\epsilon
V)^{N}\,|\Omega|^{N}\,d|\Omega|+\cdots.$

Similarly, the numerator of Eq. (55)

$\displaystyle\langle x_{0}\,e^{F}\cos\theta\rangle=\int
x_{0}\,e^{F}\,W(|\Omega|)\,\frac{1}{2^{N}N!}\,(\epsilon
V)^{N}|\Omega|^{N}\,d|\Omega|+\cdots.$ (60)

Thus, the derivative of $Z_{C}$ is obtained by

$\displaystyle-\frac{1}{V}\frac{\partial\ln
Z_{C}(T,V\rho)}{\partial\rho}\approx\frac{\int
x_{0}\,e^{F}\,W(|\Omega|)\,\frac{1}{2^{N}N!}\,(\epsilon
V)^{N}|\Omega|^{N}\,d|\Omega|}{\int
e^{F}\,W(|\Omega|)\,\frac{1}{2^{N}N!}\,(\epsilon
V)^{N}|\Omega|^{N}\,d|\Omega|}=\frac{\int x_{0}\,e^{-U_{\rm
eff}}\,d|\Omega|}{\int e^{-U_{\rm eff}}\,d|\Omega|},$ (61)

where the effective potential $U_{\rm eff}(|\Omega|)$ is defined as

$\displaystyle U_{\rm eff}(|\Omega|)=-\ln W(|\Omega|)-F-N\ln|\Omega|.$ (62)

Here, we should note that $\langle e^{F}\cos\theta\rangle$ and $\langle
x_{0}\,e^{F}\cos\theta\rangle$ go to zero in the limit of $\epsilon\to 0$,
causing the sign problem. However, the factors $(\epsilon V)^{N}$ in the
numerator and denominator of Eq. (61) are canceled, and this equation becomes
calculable without considering the additional factor $\epsilon V$. We moreover
find that this quantity is approximately equal to the value of $x_{0}$ at the
$|\Omega|$ where the effective potential $U_{\rm eff}(|\Omega|)$ is the
minimum in the case that the volume is sufficiently large.

###  V.2 Numerical simulations

Figure 3: Real part of the saddle point $x_{0}$ as a function of $|\Omega|$
for $N=10$ (blue), 50 (green), 100 (red), and 150 (black).

Figure 4: $-F$ (solid line) and $-F-N\ln|\Omega|$ (dashed line) as a function
of $|\Omega|$ for $N=10$ (blue), 50 (green), 100 (red), and 150 (black).

We perform Monte Carlo simulations of $U(1)$ lattice gauge theory with the
standard Wilson action at several inverse gauge couplings $\beta=1/g^{2}$ near
the transition point $\beta_{c}$. The effective potential $U_{\rm
eff}(|\Omega|)$ and saddle point are numerically calculated, and the effective
potential $V_{\rm eff}(N)$ for the particle number $N$ are investigated. The
lattice size is $N_{s}^{3}\times N_{t}=24^{3}\times 6$. We adopt $N_{\rm f}=2$
and $\kappa=0.025$. Since the chiral limit of free Wilson fermion is
$\kappa_{c}=1/8$, the $\kappa$ we adopt is small. The real part of saddle
point $x_{0}$ for each $|\Omega|$ and $N=\rho V$ is given by Eq. (49), which
is plotted in Fig. 4. Also, $F$ of Eq. (53) is calculated from only $|\Omega|$
and $\rho$. The solid lines in Fig. 4 are $F$ and the dashed lines are
$F-N\ln|\Omega|$ for each $|\Omega|$ and $N$. The blue, green, red, and black
curves are the results of $N=10$, 50, 100, and 150, respectively.

Using a pseudo heat bath algorithm Creutz:1983ev , the configurations are
generated at thirteen $\beta$ values:
$\beta=1.000,1.004,1.006,1.008,1.009,1.0094,1.0096,1.010,1.012,1.014,1.016,1.020,$
and $1.0240$. The data are taken until there are 1,000,000 heat bath sweeps at
each $\beta$. The multipoint reweighting method Ferrenberg:1989ui ;
Iwami:2015eba  are used to combine the data generated at different $\beta$.
The statistical errors are estimated by the jackknife method with the bin size
chosen such that the errors are saturated.

Figure 5: Polyakov-loop distribution at the transition point $\beta=1.0096$ on
the $24^{3}\times 6$ lattice.

Figure 6: $-\ln W(|\Omega|)$ as a function of $|\Omega|$ for $\beta=0.992$
(black), 1.000 (red), 1.008 (green), and 1.016 (blue).

Figure 7: Effective potential $U_{\rm eff}(|\Omega|)$ defined in Eq. (62) at
$\beta=1.00$ for $N=10$ (blue), 50 (green), 100 (red), and 150 (black).

Figure 5 shows the distribution of the Polyakov loop for each configuration at
$\beta=1.0096$. The horizontal and vertical axes are the real part and
imaginary part, respectively. This figure indicates that the phase transition
at $\beta=1.0096$ is of very weak first-order transition, where two phases
coexist. The distributions near the origin are the configurations of the
confinement phase. On the other hand, the points around the circle of
$|\Omega|=0.15$ are those of the deconfinement phase.

To compute the probability distribution function of $|\Omega|$, $W(|\Omega|)$,
we approximate the delta function with a Gaussian function like

$\displaystyle\delta(x)\approx\frac{1}{\Delta\sqrt{\pi}}\exp\left[-\left(\frac{x}{\Delta}\right)^{2}\right].$
(63)

We choose the width parameter $\Delta$ by considering a balance between the
resolution of the distribution function and its statistical error. The values
of $\Delta$ we adopt are 0.0025. The result of $-\ln W(|\Omega|)$ is shown in
Fig. 7 for $\beta=0.992$ (black), 1.000 (red), 1.008 (green), and 1.016
(blue). The results of the effective potential $U_{\rm eff}(|\Omega|)$ in Eq.
(61) at $\beta=1.00$ is plotted in Fig. 7 for $N=10$ (blue), 50 (green), 100
(red), and 150 (black), as an example. Here, since a constant may be added to
the effective potential, in Figs. 7 and 7, the constants are added so that the
minimum values of $-\ln W(|\Omega|)$ and $U_{\rm eff}(|\Omega|)$ become zero.

Figure 8: Chemical potential $\mu/T$ for which the number of particles with
the maximum generation probability is $N$ at $\kappa=0.025$ for each $\beta$
in $U(1)$ lattice gauge theory with two flavors.

We then calculate the derivative of the canonical partition function. We use
the following equation to remove $W(|\Omega|)$ from Eq. (61) so that the
equation is independent of the approximation of the delta function, Eq. (63),

$\displaystyle-\frac{1}{V}\frac{\partial\ln
Z_{C}(T,V\rho)}{\partial\rho}=\frac{\int{\cal
D}U\,x_{0}\,e^{F}|\hat{\Omega}|^{N}}{\int{\cal D}U\,e^{F}|\hat{\Omega}|^{N}}$
(64)

This quantity is $\mu/T$ where the density with the maximum generation
probability is $\rho$ and is computable without the sign problem. We plot the
result as a function of $N=V\rho$ for $\beta=0.992$ (purple), 0.996 (magenta),
1.00 (red), 1.004 (blue), 1.008 (green), 1.0012 (black), and 1.016 (brown) in
Fig. 8. The fine up and down vibrations in Fig. 8 may be an artifact caused by
the fact that the effective potential in Fig. 7 is not smooth, since this
quantity is the value of $x_{0}$ at the $|\Omega|$ where $U_{\rm eff}$ is the
minimum. Ignoring the fine vibrations, when $\beta$ is in the deconfinement
phase at zero density, the chemical potential is monotonically increasing.
However, in the case of a confined phase at zero density, as the density
increases, the chemical potential drops once and increases again. This means
that there are multiple $N$s for a specific $\mu/T$. This is a typical sign
when crossing the first-order phase transition.

###  V.3 Higher-order terms of the hopping parameter expansion

As seen in the previous section, our method to solve the problem of the center
symmetry and the sign problem works well in $U(1)$ gauge theory. In the above
calculation, the quark determinant has been evaluated by approximating only
the leading-order term of the hopping parameter expansion. We extend this
method to include higher-order terms of the expansion.

We classify each term of the hopping parameter expansion by winding number and
denote $D(z)$ in Eq. (38) as follows,

$\displaystyle D(z)=G_{0}+\sum_{n=1}^{n_{\rm
max}}\left[e^{nz}G_{n}+e^{-nz}G_{n}^{*}\right]$ (65)

with $n_{\rm max}=\infty$, where $G_{n}$ is the sum of terms whose winding
number is $n$ in the positive direction, and $G_{n}^{*}$ is those of winding
number $n$ in the negative direction. The first and second derivatives of
$D(z)$ with respect to $z$ are

$\displaystyle\frac{\partial D(z)}{\partial z}=\sum_{n=1}^{n_{\rm
max}}\left[ne^{nz}G_{n}-ne^{-nz}G_{n}^{*}\right],\hskip
14.22636pt\frac{\partial^{2}D(z)}{\partial z^{2}}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{n_{\rm
max}}\left[n^{2}e^{nz}G_{n}+n^{2}e^{-nz}G_{n}^{*}\right].$ (66)

The saddle point $z_{0}=x_{0}+iy_{0}$ is defined as

$\displaystyle\frac{\partial D}{\partial z}(z_{0})=\sum_{n=1}^{n_{\rm
max}}\left[ne^{nz_{0}}G_{n}-ne^{-nz_{0}}G_{n}^{*}\right]=\rho.$ (67)

For the case that $D(z)$ can be approximated by $n_{\rm max}=1$, i.e.
$D(z)=G_{0}+e^{z}G_{1}+e^{-z}G_{1}^{*}$. The analysis is almost the same as
the leading-order analysis. Because $\rho$ is a real number, the complex phase
of $(e^{iy_{0}}G_{1}-e^{-iy_{0}}G_{1}^{*})$ is zero. Then, the imaginary part
of the saddle point $y_{0}$ is determined as

$\displaystyle\tan y_{0}=-\frac{{\rm Im}G_{1}}{{\rm Re}G_{1}}=-{\rm Arg}\
G_{1}.$ (68)

Then, the real part $x_{0}$ is given by

$\displaystyle e^{x_{0}}|G_{1}|-e^{-x_{0}}|G_{1}|=\rho.$ (69)

The real part of the saddle point $x_{0}$ is

$\displaystyle x_{0}=\mathrm{arcsinh}\left(\frac{\rho}{2|G_{1}|}\right).$ (70)

$D(z_{0})$ and $\partial^{2}D/\partial z^{2}(z_{0})$ at the saddle point are
real numbers,

$\displaystyle D(z_{0})=F_{0}+(e^{x_{0}}+e^{-x_{0}})|G_{1}|,\hskip
14.22636pt\frac{\partial^{2}D}{\partial
z^{2}}(z_{0})=(e^{x_{0}}+e^{-x_{0}})|G_{1}|.$ (71)

Substituting $(x_{0},y_{0})$, $D(z_{0})$ and $\partial^{2}D/\partial
z^{2}(z_{0})$ into Eqs. (53) and (54). $-\partial\ln
Z_{C}(T,V\rho)/\partial\rho$ can be computed as well as the leading-order
calculation.

The convergence of the hopping parameter expansion is investigated in Ref.
Wakabayashi:2021eye . In the hopping parameter expansion, the nonzero
contribution of $G_{n}$ appears from $O(\kappa^{nN_{t}})$ and $G_{n}$ with
$n\geq 2$ is much smaller than $G_{1}$ in relatively low-order terms of
$\kappa$. Thus, $D(z)$ is well approximated by $n_{\rm max}=1$ for relatively
small $\kappa$. In such cases, the sign problem does not appear at all.
However, as the quark mass decreases, the contributions from higher-order
terms of $\kappa$-expansion increase. Therefore, it is necessary to consider
$G_{n}$ with $n\geq 2$. The complex phase for $n_{\rm max}\geq 2$ is not as
easy to remove as when $n_{\rm max}=1$. Then, the complex phase of $D(z_{0})$
and $\partial^{2}D/\partial z^{2}(z_{0})$ must be taken into account.

##  VI Application to $SU(3)$ lattice gauge theory

In the discussion of avoiding sign problems in $U(1)$ gauge theory, the $U(1)$
center symmetry is essentially important, while the center symmetry is $Z_{3}$
for $SU(3)$ gauge theory. Here, we discuss the probability distribution of the
Polyakov loop on the complex plane in $SU(3)$ gauge theory. To clarify this
argument, we should consider the Polyakov loop without spatial averaging
(local Polyakov loop) and their spatial average $\Omega$ (averaged Polyakov
loop) separately.

###  VI.1 Distribution of local Polyakov loop

First, we discuss the distribution of the Polyakov loop at one spatial point.
We define an element of $SU(3)$ group as $U=e^{iH}$. $H$ is a traceless
Hermitian matrix, i.e. $\mathrm{tr}H=0$ and $H=H^{{\dagger}}$. We diagonalize
$H$ as $H=V\Lambda V^{{\dagger}}$, where $V$ is a unitary matrix and $\Lambda$
is a diagonal matrix with the diagonal elements
$(\lambda_{1},\lambda_{2},\lambda_{3})$. Because $H=H^{{\dagger}}$,
$\lambda_{1},\lambda_{2}$ and $\lambda_{3}$ are real numbers, and
$\mathrm{tr}H=\mathrm{tr}\Lambda=\lambda_{1}+\lambda_{2}+\lambda_{3}=0$. Since
the Polyakov loop is the trace of a $SU(3)$ matrix divided by 3,
$\mathrm{tr}U/3,$ the trace can be rewritten as

$\displaystyle\mathrm{tr}U$ $\displaystyle=$
$\displaystyle\mathrm{tr}\left(e^{iH}\right)=\mathrm{tr}\sum_{n=1}^{\infty}\frac{(iH)^{n}}{n!}=\sum_{n=1}^{\infty}\frac{i^{n}}{n!}\mathrm{tr}\left(H^{n}\right)=\sum_{n=1}^{\infty}\frac{i^{n}}{n!}\mathrm{tr}\left(\Lambda^{n}\right)$
(72) $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{i^{n}}{n!}(\lambda_{1}^{n}+\lambda_{2}^{n}+\lambda_{3}^{n})=e^{i\lambda_{1}}+e^{i\lambda_{2}}+e^{i\lambda_{3}}.$

The complex number that can be the value of $\mathrm{tr}U$ is in the triangle
drawn by the red lines in Fig. 9. The square of the absolute value of
$\mathrm{tr}U$ is given as

$\displaystyle|\mathrm{tr}U|^{2}$ $\displaystyle=$
$\displaystyle(e^{i\lambda_{1}}+e^{i\lambda_{2}}+e^{i\lambda_{3}})(e^{-i\lambda_{1}}+e^{-i\lambda_{2}}+e^{-i\lambda_{3}})$
(73) $\displaystyle=$ $\displaystyle
3+e^{i(\lambda_{1}-\lambda_{2})}+e^{i(\lambda_{1}-\lambda_{3})}+e^{i(\lambda_{2}-\lambda_{1})}+e^{i(\lambda_{2}-\lambda_{3})}+e^{i(\lambda_{3}-\lambda_{1})}+e^{i(\lambda_{3}-\lambda_{2})}$
$\displaystyle=$ $\displaystyle
3+2\cos(\lambda_{1}-\lambda_{2})+2\cos(\lambda_{2}-\lambda_{3})+2\cos(\lambda_{3}-\lambda_{1}).$

The condition for maximizing $|\mathrm{tr}U|^{2}$ is
$\lambda_{1}=\lambda_{2}=\lambda_{3}$ in $\mathrm{mod}\ 2\pi$. Because
$\lambda_{1}+\lambda_{2}+\lambda_{3}=0$,
$\lambda_{1}=\lambda_{2}=\lambda_{3}=0$, $2\pi/3$ or $4\pi/3$ in
$\mathrm{mod}\ 2\pi$. Then, $|\mathrm{tr}U|=3$, and the Polyakov loop is
$\mathrm{tr}U/3=(1,0),(-1/2,\sqrt{3}/2)$, or $(-1/2,-\sqrt{3}/2)$. The
condition for the boundary is $\lambda_{1}=\lambda_{2}\neq\lambda_{3}$,
$\lambda_{2}=\lambda_{3}\neq\lambda_{1}$ or
$\lambda_{3}=\lambda_{1}\neq\lambda_{2}$. In the case of
$\lambda_{1}=\lambda_{2}\neq\lambda_{3}$, for example,
$|\mathrm{tr}U|^{2}=5+4\cos x$ as a function of $x=\lambda_{3}-\lambda_{1}$.
Then, $\mathrm{tr}U/3=(2e^{-ix/3}+e^{2ix/3})/3$, which is the red curve in
Fig. 9.

Figure 9: The distribution of the local Polyakov loop in a complex plane on
one configuration of $SU(3)$ gauge theory generated at $\beta=5.60$ (left) and
$\beta=6.00$ (right).

We perform simulations of $SU(3)$ gauge theory (quenched QCD) on a lattice
with the size $24^{3}\times 4$. The phase transition point is
$\beta_{c}=5.69138(3)$ on the $24^{3}\times 4$ lattice Saito:2011fs . Focusing
on one configuration after sufficient thermalization, we investigate the
distribution of the local Polyakov loop. In the left and right panels of Fig.
9, we plot the local Polyakov loop at each point on a configuration in the
confinement phase $\beta=5.6$ and the deconfinement phase $\beta=6.0$,
respectively. The distributions are in the triangle drown by red curves. In
the confinement phase, the distribution of the local Polyakov loops at each
point is $Z_{3}$ symmetric. On the other hand, the $Z_{3}$ symmetry in the
distribution of the local Polyakov loop is broken in the deconfinement phase.

###  VI.2 Distribution of averaged Polyakov loop

Next, we discuss the distribution of the averaged Polyakov loop. We map the
real and imaginary parts of the local Polyakov loop to two-dimensional plane
of $\vec{x}=(x,y)=(r\cos\theta,r\sin\theta)$. We calculate the variance in a
certain direction. Let $\vec{e}$ be the unit vector in the direction
perpendicular to that direction, and let $\eta$ be the angle between $\vec{e}$
and the $x$-axis. The distance from the straight line drawn in the $\vec{e}$
direction from the origin is $|\vec{x}\times\vec{e}|=r\sin(\theta-\eta)$. The
variance is $\int|\vec{x}\times\vec{e}|^{2}P(\vec{x})\ d\vec{x}$, where
$P(\vec{x})$ is the probability that the complex value of the local Polyakov
loop will be $\vec{x}$. Here, we assume that the probability distribution is
$Z_{N}$ symmetric. Then, the generation probabilities of $Z_{N}$ symmetric
points: $\vec{x}_{i}$ with $i=1,2,\cdots,N$ satisfying
$\vec{x}_{1}+\vec{x}_{2}+\cdots+\vec{x}_{N}=0$ and
$|\vec{x}_{1}|=|\vec{x}_{2}|=\cdots=|\vec{x}_{N}|=r$, are equal, i.e.
$P(\vec{x}_{1})=P(\vec{x}_{2})=\cdots=P(\vec{x}_{N})$. The variance is

$\displaystyle\int|\vec{x}\times\vec{e}|^{2}P(\vec{x})\
d\vec{x}=\frac{1}{N}\sum_{i=1}^{N}\int|\vec{x}_{i}\times\vec{e}|^{2}P(\vec{x}_{i})\
d\vec{x}_{i}=\int\frac{1}{N}\sum_{i=1}^{N}|\vec{x}_{i}\times\vec{e}|^{2}P(\vec{x}_{1})\
d\vec{x}_{1}.$ (74)

We calculate the sum for symmetric points before integrating over $\vec{x}$,

$\displaystyle\chi^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}|\vec{x}_{i}\times\vec{e}|^{2}=\frac{1}{N}\sum_{i=1}^{N}r^{2}\sin^{2}(\theta_{i}-\eta)=\frac{r^{2}}{N}\sum_{i=1}^{N}\frac{1-\cos
2(\theta_{i}-\eta)}{2}$ (75) $\displaystyle=$
$\displaystyle\frac{r^{2}}{N}\left(\frac{N}{2}-\frac{1}{2}\sum_{i=1}^{N}\cos
2(\theta_{i}-\eta)\right)=r^{2}\left(\frac{1}{2}-\frac{1}{2N}\mathrm{Re}\left[e^{-2i\eta}\sum_{i=1}^{N}e^{2i\theta_{i}}\right]\right),$

where $\vec{x}_{i}=(x_{i},y_{i})=(r\cos\theta_{i},r\sin\theta_{i})$. Because
$\theta_{n}$ is given by $\theta_{n}=\theta_{1}+2\pi(n-1)/N$,

$\displaystyle\sum_{n=1}^{N}e^{2i\theta_{n}}=\sum_{n=1}^{N}e^{2i(\theta_{1}+2\pi(n-1)/N)}=e^{2i\theta_{1}}\sum_{n=1}^{N}e^{4i\pi(n-1)/N}=0,$
(76)

except for $N=2$. Hence, $\chi^{2}=r^{2}/2=|\vec{x}_{1}|^{2}/2$ for $N\geq 3$,
and

$\displaystyle\int|\vec{x}\times\vec{e}|^{2}P(\vec{x})\
d\vec{x}=\int\frac{|\vec{x}|^{2}}{2}P(\vec{x})\ d\vec{x}.$ (77)

The result of the variance does not depend on $\eta$. This means that the
variances in the real axis direction and the imaginary axis direction are the
same when the probability distribution is $Z_{N}$ symmetric.

For the case that the distribution of the local Polyakov loop has $Z_{N}$
symmetry with $N\geq 3$, the variances of the real and imaginary parts are the
same. When the volume is large enough, the distributions of the real and
imaginary parts of the averaged Polyakov loop are Gaussian by the central
limit theorem. The width of the Gaussian distribution is the same for the real
and imaginary parts. Then, the probability distribution of the averaged
Polyakov loop $\Omega=|\Omega|e^{i\phi}$ is

$\displaystyle W(|\Omega|,\phi)\approx
Ce^{-\alpha(\mathrm{Re}\Omega)^{2}}e^{-\alpha(\mathrm{Im}\Omega)^{2}}=Ce^{-\alpha|\Omega|^{2}},$
(78)

where the parameter $\alpha$ is inversely proportional to the variance of the
local Polyakov loop distribution and $C$ is the normalization constant. The
distribution does not depend on the complex phase $\phi$ and has $U(1)$
symmetry.

To summarize the above discussion, the distribution of the local Polyakov
loops is $Z_{3}$ symmetric in the confinement phase. Then, the probability
distribution of the spatially averaged Polyakov loop $\Omega$ is $U(1)$
symmetric in the volume infinity limit, which is not $Z_{3}$ symmetric. Since
the distribution of $\Omega$ is $U(1)$ symmetric, the method of avoiding the
sign problem used in $U(1)$ gauge theory can be applied. On the other hand, in
the deconfinement phase of $SU(3)$ gauge theory, $Z_{3}$ symmetry is
spontaneously broken as seen in the right panel of Fig. 9. Then, the
probability distribution of the averaged Polyakov loop is not $U(1)$
symmetric. Thus, the method to avoid the sign problem in $U(1)$ gauge theory
cannot be applied. However, the sign problem is not serious in the
deconfinement phase. The complex phase distribution can be well approximated
by the Gaussian distribution. In such a case, the sign problem may be avoided
by the method used in Refs. Ejiri:2007ga ; Ejiri:2009hq ; Nakagawa:2011eu ;
Nakagawa:2012yx ; Saito:2013vja .

##  VII Conclusions

We studied the probability distribution function of particle density. The
probability distribution function is obtained by calculating the canonical
partition function fixing the number of particles from the grand partition
function. However, if the system has the center symmetry on a finite lattice,
the canonical partition function is zero when the number of particles is not a
multiple of three for $SU(3)$ gauge theory. For $U(1)$ gauge theory, the
canonical partition function is zero except when the particle number is zero.
Then, the probability distribution function is zero for these cases. This
situation is natural in the confined phase, but is unacceptable in the
deconfinement phase because there should be states of various particle
numbers.

This problem is essentially the same as the problem that the expectation value
of the Polyakov loop is always zero when calculating with finite volume, due
to the center symmetry. To solve this problem, it is necessary to add an
infinitesimal external field to break the symmetry and take the limit of
infinite volume. Moreover, in the case of $U(1)$ gauge theory, the sign
problem can be solved using the $U(1)$ center symmetry at the same time.

We performed numerical simulations of $U(1)$ lattice gauge theory near the
deconfinement phase transition point. When the dynamical fermions are heavy,
we actually demonstrated that the calculation of the probability distribution
function at finite density is possible using the method proposed in this
study. We calculated the derivative of the canonical partition function using
a saddle point approximation Ejiri:2008xt , and found that our method to avoid
the sign problem works well. From the canonical partition function, we
calculated $\mu/T$ as a function of density $\rho$. Then, the nature of the
phase transition can be investigated.

The application of this method to QCD ($SU(3)$ gauge theory) was discussed. If
the distribution of the local Polyakov loops is $Z_{3}$ symmetric in the
confinement phase, the probability distribution of the spatially averaged
Polyakov loop $\Omega$ is $U(1)$ symmetric when the spatial volume is
sufficiently large. Then, the method of avoiding the sign problem and solving
the problem of the center symmetry used in $U(1)$ gauge theory can be applied.
On the other hand, in the deconfinement phase, the probability distribution of
the averaged Polyakov loop is not $U(1)$ symmetric. Thus, the method to avoid
the sign problem in $U(1)$ gauge theory cannot be applied. However, the sign
problem is not very serious in the deconfinement phase.

### Acknowledgments

The author thanks the members of the WHOT-QCD Collaboration for useful
discussions. This work was supported by JSPS KAKENHI Grants No. JP21K03550,
No. JP20H01903, and No. JP19H05146, and the Uchida Energy Science Promotion
Foundation. The author also thanks the Yukawa Institute for Theoretical
Physics at Kyoto University for the workshops YITP-W-19-09.

## References

  * (1) S. Ejiri, F. Karsch and K. Redlich, “Hadronic fluctuations at the QCD phase transition,” Phys. Lett. B 633, 275-282 (2006).
  * (2) M. A. Stephanov, “Non-Gaussian fluctuations near the QCD critical point,” Phys. Rev. Lett. 102, 032301 (2009).
  * (3) M. Asakawa, S. Ejiri and M. Kitazawa, “Third moments of conserved charges as probes of QCD phase structure,” Phys. Rev. Lett. 103, 262301 (2009)
  * (4) M. Cheng, P. Hegde, C. Jung, F. Karsch, O. Kaczmarek, E. Laermann, R. D. Mawhinney, C. Miao, P. Petreczky and C. Schmidt, W. Soeldner, “Baryon Number, Strangeness and Electric Charge Fluctuations in QCD at High Temperature,” Phys. Rev. D 79, 074505 (2009).
  * (5) C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and K. Redlich, “Thermodynamics of two flavor QCD to sixth order in quark chemical potential,” Phys. Rev. D 71, 054508 (2005).
  * (6) C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and C. Schmidt, “The Equation of state for two flavor QCD at nonzero chemical potential,” Phys. Rev. D 68, 014507 (2003).
  * (7) S. Ejiri, Y. Maezawa, N. Ukita, S. Aoki, T. Hatsuda, N. Ishii, K. Kanaya, T. Umeda [WHOT-QCD], “Equation of State and Heavy-Quark Free Energy at Finite Temperature and Density in Two Flavor Lattice QCD with Wilson Quark Action,” Phys. Rev. D 82, 014508 (2010).
  * (8) M. D’Elia, G. Gagliardi and F. Sanfilippo, “Higher order quark number fluctuations via imaginary chemical potentials in $N_{f}=2+1$ QCD,” Phys. Rev. D 95, 094503 (2017).
  * (9) A. Bazavov, H. T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa, S. Mukherjee, H. Ohno and P. Petreczky, et al. “The QCD Equation of State to $\mathcal{O}(\mu_{B}^{6})$ from Lattice QCD,” Phys. Rev. D 95, 054504 (2017).
  * (10) S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K. Szabo, A. Pasztor, I. Portillo and C. Ratti, “Higher order fluctuations and correlations of conserved charges from lattice QCD,” JHEP 10, 205 (2018).
  * (11) A. Bazavov, D. Bollweg, H. T. Ding, P. Enns, J. Goswami, P. Hegde, O. Kaczmarek, F. Karsch, R. Larsen and S. Mukherjee, et al. “Skewness, kurtosis, and the fifth and sixth order cumulants of net baryon-number distributions from lattice QCD confront high-statistics STAR data,” Phys. Rev. D 101, 074502 (2020).
  * (12) R. Bellwied, S. Borsanyi, Z. Fodor, J. N. Guenther, S. D. Katz, P. Parotto, A. Pasztor, D. Pesznyak, C. Ratti and K. K. Szabo, “Corrections to the hadron resonance gas from lattice QCD and their effect on fluctuation-ratios at finite density,” Phys. Rev. D 104, 094508 (2021).
  * (13) A. Hasenfratz and D. Toussaint, “Canonical ensembles and nonzero density quantum chromodynamics,” Nucl. Phys. B 371, 539-549 (1992).
  * (14) P. de Forcrand and S. Kratochvila, “Finite density QCD with a canonical approach,” Nucl. Phys. B Proc. Suppl. 153, 62-67 (2006).
  * (15) S. Kratochvila and P. de Forcrand, “The Canonical approach to finite density QCD,” PoS LAT2005, 167 (2005).
  * (16) S. Ejiri, “Canonical partition function and finite density phase transition in lattice QCD,” Phys. Rev. D 78, 074507 (2008).
  * (17) A. Alexandru, M. Faber, I. Horvath and K. F. Liu, “Lattice QCD at finite density via a new canonical approach,” Phys. Rev. D 72, 114513 (2005).
  * (18) A. Li, A. Alexandru, K. F. Liu and X. Meng, “Finite density phase transition of QCD with $N_{f}=4$ and $N_{f}=2$ using canonical ensemble method,” Phys. Rev. D 82, 054502 (2010).
  * (19) A. Li, A. Alexandru and K. F. Liu, “Critical point of $N_{f}=3$ QCD from lattice simulations in the canonical ensemble,” Phys. Rev. D 84, 071503 (2011).
  * (20) J. Danzer and C. Gattringer, “Properties of canonical determinants and a test of fugacity expansion for finite density lattice QCD with Wilson fermions,” Phys. Rev. D 86, 014502 (2012).
  * (21) R. Fukuda, A. Nakamura and S. Oka, “Canonical approach to finite density QCD with multiple precision computation,” Phys. Rev. D 93, 094508 (2016).
  * (22) A. Nakamura, S. Oka and Y. Taniguchi, “QCD phase transition at real chemical potential with canonical approach,” JHEP 02, 054 (2016).
  * (23) V. G. Bornyakov, D. L. Boyda, V. A. Goy, A. V. Molochkov, A. Nakamura, A. A. Nikolaev and V. I. Zakharov, “New approach to canonical partition functions computation in $N_{f}=2$ lattice QCD at finite baryon density,” Phys. Rev. D 95, 094506 (2017)
  * (24) M. Wakayama, V. G. Bornyakov, D. L. Boyda, V. A. Goy, H. Iida, A. V. Molochkov, A. Nakamura and V. I. Zakharov, “Lee-Yang zeros in lattice QCD for searching phase transition points,” Phys. Lett. B 793, 227-233 (2019)
  * (25) I. Montvay G. Munster, “Quantum Fields on a Lattice,” Cambridge University Press (1994)
  * (26) H. J. Rothe, “Lattice Gauge Theories: An Introduction”, World Scientific (1992)
  * (27) H. Saito, S. Ejiri, S. Aoki, T. Hatsuda, K. Kanaya, Y. Maezawa, H. Ohno, T. Umeda [WHOT-QCD], “Phase structure of finite temperature QCD in the heavy quark region,” Phys. Rev. D 84, 054502 (2011).
  * (28) H. Saito, S. Ejiri, S. Aoki, K. Kanaya, Y. Nakagawa, H. Ohno, K. Okuno and T. Umeda [WHOT-QCD], “Histograms in heavy-quark QCD at finite temperature and density,” Phys. Rev. D 89, 034507 (2014).
  * (29) N. Wakabayashi, S. Ejiri, K. Kanaya and M. Kitazawa, “Scope and convergence of the hopping parameter expansion in finite temperature QCD with heavy quarks around the critical point,” PTEP 2022, 033B05 (2022).
  * (30) S. Ejiri, S. Itagaki, R. Iwami, K. Kanaya, M. Kitazawa, A. Kiyohara, M. Shirogane, T. Umeda [WHOT-QCD], “End point of the first-order phase transition of QCD in the heavy quark region by reweighting from quenched QCD,” Phys. Rev. D 101, 054505 (2020).
  * (31) A. Kiyohara, M. Kitazawa, S. Ejiri and K. Kanaya, “Finite-size scaling around the critical point in the heavy quark region of QCD,” Phys. Rev. D 104, 114509 (2021).
  * (32) A. Roberge and N. Weiss, “Gauge Theories With Imaginary Chemical Potential and the Phases of QCD,” Nucl. Phys. B 275, 734-745 (1986).
  * (33) M. Creutz, L. Jacobs and C. Rebbi, “Monte Carlo Computations in Lattice Gauge Theories,” Phys. Rept. 95, 201-282 (1983).
  * (34) A. M. Ferrenberg and R. H. Swendsen, “Optimized Monte Carlo analysis,” Phys. Rev. Lett. 63, 1195-1198 (1989).
  * (35) R. Iwami, S. Ejiri, K. Kanaya, Y. Nakagawa, D. Yamamoto and T. Umeda, “Multipoint reweighting method and its applications to lattice QCD,” Phys. Rev. D 92, 094507 (2015).
  * (36) S. Ejiri, “On the existence of the critical point in finite density lattice QCD,” Phys. Rev. D 77, 014508 (2008).
  * (37) Y. Nakagawa, S. Ejiri, S. Aoki, K. Kanaya, H. Ohno, H. Saito, T. Hatsuda, T. Umeda [WHOT-QCD], “Histogram method in finite density QCD with phase quenched simulations,” PoS LATTICE2011, 208 (2011).
  * (38) Y. Nakagawa, S. Ejiri, S. Aoki, K. Kanaya, H. Saito, H. Ohno, T. Hatsuda, T. Umeda [WHOT-QCD], “Phase structure of finite density QCD with a histogram method,” PoS LATTICE2012, 092 (2012).