# Finite-size scaling around the critical point in the heavy quark region of
QCD

Atsushi Kiyohara  Department of Physics, Osaka University, Toyonaka, Osaka
560-0043, Japan  Masakiyo Kitazawa  Department of Physics, Osaka University,
Toyonaka, Osaka 560-0043, Japan  J-PARC Branch, KEK Theory Center, Institute
of Particle and Nuclear Studies, KEK, 203-1, Shirakata, Tokai, Ibaraki
319-1106, Japan  Shinji Ejiri  Department of Physics, Niigata University,
Niigata 950-2181, Japan  Kazuyuki Kanaya  Tomonaga Center for the History of
the Universe, University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan

###### Abstract

Finite-size scaling is investigated in detail around the critical point in the
heavy-quark region of nonzero temperature QCD. Numerical simulations are
performed with large spatial volumes up to the aspect ratio $N_{s}/N_{t}=12$
at a fixed lattice spacing with $N_{t}=4$. We show that the Binder cumulant
and the distribution function of the Polyakov loop follow the finite-size
scaling in the $Z(2)$ universality class for large spatial volumes with
$N_{s}/N_{t}\geq 9$, while, for $N_{s}/N_{t}\leq 8$, the Binder cumulant
becomes inconsistent with the $Z(2)$ scaling. To realize the large-volume
simulations in the heavy-quark region, we adopt the hopping tails expansion
for the quark determinant: We generate gauge configurations using the leading
order action including the Polyakov loop term for $N_{t}=4$, and incorporate
the next-to-leading order effects in the measurements by the multipoint
reweighting method. We find that the use of the leading-order configurations
is crucially effective in suppressing the overlapping problem in the
reweighting and thus reducing the statistical errors.

††preprint: J-PARC-TH-0246,UTHEP-758

##  I Introduction

One of the interesting features of the medium described by quantum
chromodynamics (QCD) is the existence of phase transitions of various orders.
While the finite-temperature QCD transition is an analytic crossover at zero
quark chemical potential $\mu_{q}$ [1, 2], this phase transition is expected
to become of first order in dense medium with large $\mu_{q}$ [3]. The end
point of the first-order transition is called the critical point (CP) at which
the transition is of second order. The singularity in thermodynamic
observables associated with the second order CP are believed to be useful in
detecting the CP in heavy-ion collision experiments [4, 5]. Accordingly,
researches called the beam-energy scan are actively performed in experimental
facilities all over the world to search for the critical fluctuations around
the CP [5, 6, 7].

The order of the finite temperature QCD transition changes also with variation
of quark masses [8]. For the $2+1$-flavor QCD, it is known that the crossover
at the physical quark masses becomes of first order both in the light- and
heavy-quark limits; the phase diagram representing this feature is known as
the Columbia plot [9, 10]. Revealing the nature of the phase transitions with
the variation of quark masses is an important subject of QCD at nonzero
temperature since it provides us with various insights into the transition at
the physical quark masses.

Pinning down the boundaries of the first-order transitions in $2+1$-flavor QCD
in the light [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] and heavy
[24, 25, 26, 27, 28, 29, 30, 31] quark regions is a longstanding subject in
lattice QCD simulations. It, however, has been found that the location of the
boundaries are strongly dependent on the lattice cutoff of the simulations
[22, 23, 31], and their quantitative determination in the continuum limit has
not been established yet.

One of the difficulties in these analyses is that observables near the CP are
strongly dependent on the spatial volume of the system. The spatial volume
dependence is in part described by the finite-size scaling (FSS) [32].
However, the FSS is applicable only for describing the singular part of
thermodynamic quantities that dominates over the non-singular part only in the
vicinity of the CP for sufficiently large spatial volumes. When the spatial
volume is not large enough the FSS of observables is violated due to the
contributions of the non-singular part and this makes their analysis based on
the FSS problematic. In fact, although the CP of QCD is believed to belong to
the three-dimensional $Z(2)$ universality class [8], a clear FSS in this
universality class has not been observed in latest numerical studies around
the CP in the light [23] and heavy [31] quark regions. These results suggest
the necessity to perform numerical analyses with yet larger spatial volumes.

In the present study, we focus on the CP in the heavy quark region and study
the behavior of observables by numerical simulations with large spatial
volumes corresponding to the aspect ratios up to $N_{s}/N_{t}=12$. To carry
out analysis on the large spatial volumes with high precision, we fix the
temporal lattice extent to be $N_{t}=4$ in this study.

We also employ the hopping parameter expansion (HPE) to deal with the quark
determinant. In this study, we generate gauge configurations using the leading
order (LO) action of the HPE including the Polyakov loop term for $N_{t}=4$,
and then incorporate the next-to-leading order (NLO) effects by a multipoint
reweighting method [33]. In our previous study, we generated the
configurations in quenched QCD (zero-th order of the HPE) and reweighted them
to incorporate LO and NLO effects [27, 29, 30]. We show that the use of the LO
action to generate configurations is quite effective in suppressing the
overlapping problem of the reweighting method. This is essential to carry out
simulations with large system volumes as studied in the present paper. We also
verify the convergence of the HPE by comparing the LO and NLO results.

We perform the Binder cumulant analysis of the Polyakov loop around the CP and
find that the numerical results for $N_{s}/N_{t}\geq 9$ follow well the FSS in
the $Z(2)$ universality class. On the other hand, inclusion of the data at
$N_{s}/N_{t}\leq 8$ in the analysis gives rise to a statistically-significant
deviation from the scaling behavior. We further investigate the scaling
behavior of the distribution function of the Polyakov loop. We find that the
distribution function follows the FSS in the $Z(2)$ universality class for
$N_{s}/N_{t}\geq 8$. From the deviation pattern of the distribution function
for $LT=6$ from the $Z(2)$ FSS, we discuss that the violation of the scaling
behavior in the Binder cumulant is caused by the deviation in the tails of the
distribution for small $LT$.

This paper is organized as follows. In the next section we give a brief review
on the FSS. We then explain the setup of our lattice simulation and analyses
using the HPE in Sec. III. In Sec. IV, we determine the transition line and
perform the Binder cumulant analysis to determine the location of the CP and
the critical exponent. In Sec. V, we investigate the FSS of the distribution
function of the Polyakov loop. The last section, Sec. VI, is devoted to a
summary. In Appendix A, we give definition of cumulants. In Appendix B, the
HPE of the quark determinant is calculated up to the NLO. In Appendix C, the
convergence of the HPE is examined by comparing the Binder cumulants at the LO
and the NLO.

##  II Finite-size scaling

Figure 1:  Phase diagrams of (a) QCD in heavy-quark region and (b) three-
dimensional Ising model.

Let us first give a brief review on the FSS and its application to the CP in
the heavy-quark region of QCD.

The heavy-quark limit of QCD corresponds to the $SU(3)$ Yang-Mills theory
(quenched QCD). This theory has a first-order deconfinement phase transition
at nonzero temperature $T$. When the quark mass $m_{q}$ is finite (throughout
this section we assume that the quark masses are degenerate), with increasing
$1/m_{q}$, this first-order transition becomes weaker and eventually
terminates at the CP, as schematically shown in the phase diagram on the
$(T,1/m_{q})$ plane in Fig. 1 (a). This CP, as well as that in the light quark
region, is believed to belong to the $Z(2)$ universality class, i.e. the
universality class of the three-dimensional Ising model [8].

Near the CP of the three-dimensional Ising model, the relevant scaling
parameters are the reduced temperature $t$ and external magnetic field $h$;
extensive variables conjugate to these parameters are the energy and the
magnetization, respectively. As shown in Fig. 1 (b), the CP is located at
$(t,h)=(0,0)$ and the first-order transition exists on the $t$ axis for $t<0$.
The singular part of thermodynamic quantities near the CP is described by the
scaling function of $t$ and $h$. According to the universality, the singular
part of thermodynamic quantities near the CP of heavy-quark QCD is described
by the same scaling function, where the scaling parameters $t$ and $h$ are
encoded into the $(T,1/m_{q})$ plane as schematically shown in Fig. 1 (a); the
$t$ axis is parallel to the first-order line at the CP while the direction of
the $h$ axis is not constrained from the universality.

According to the FSS argument [32] the singular part of the dimension-less
free energy $F(t,h,L^{-1})$ around the CP obtained at a finite volume
$V=L^{3}$ has a scaling

$\displaystyle F(t,h,L^{-1})=F(tb^{y_{t}},hb^{y_{h}},L^{-1}b),$ (1)

for arbitrary scale factor $b$. The values of the exponents $y_{t}$ and
$y_{h}$ are specific for each universality class. In the $Z(2)$ universality
class these parameters are numerically obtained as [32]

$\displaystyle y_{t}=0.110,\quad y_{h}=-1.237.$ (2)

By setting $b=L$ one has

$\displaystyle
F(t,h,L^{-1})=F(tL^{y_{t}},hL^{y_{h}},1)\equiv\tilde{F}(tL^{y_{t}},hL^{y_{h}}).$
(3)

Derivatives of $F(t,h,L^{-1})$ with respect to $t$ and $h$ define the
cumulants of the corresponding extensive variables. For example, cumulants of
the magnetization $M$ are given by

$\displaystyle\langle M(t,h,L^{-1})^{n}\rangle_{\rm
c}=\partial_{h}^{n}F(t,h,L^{-1}),$ (4)

with $\partial_{h}=\partial/\partial h$. From Eq. (3), $L$ dependence of the
cumulants $\langle M^{n}\rangle_{\rm c}$ near the CP are written as

$\displaystyle\langle M(t,h,L^{-1})^{n}\rangle_{\rm c}$
$\displaystyle=\partial_{h}^{n}F(t,h,L^{-1})$
$\displaystyle=L^{ny_{h}}\partial_{h}^{n}\tilde{F}(tL^{y_{t}},hL^{y_{h}}).$
(5)

In Ref. [34], it is suggested that the so-called (fourth-order) Binder
cumulant

$\displaystyle B_{4}(t,h,L^{-1})$ $\displaystyle=\frac{\langle
M(t,h,L^{-1})^{4}\rangle_{\rm c}}{(\langle M(t,h,L^{-1})^{2}\rangle_{\rm
c})^{2}}+3$ (6)

plays a useful role to determine the location of the CP from numerical results
obtained at finite $L$. When the distribution of $M$ obeys the Gauss
distribution or the distribution composed of two delta functions with an equal
weight, we have

$\displaystyle B_{4}=\begin{cases}3&\mathrm{Gauss~{}distribution},\\\
1&\mathrm{two~{}delta~{}functions},\end{cases}$ (7)

respectively. Since the distribution of $M$ approaches these functions in the
$L\to\infty$ limit on the crossover and first-order lines at $h=0$,
respectively, $B_{4}$ should approach Eq. (7) on these lines in the
$L\to\infty$ limit. Moreover, from Eq. (4), $B_{4}$ at $h=0$ behaves as a
function of $t$ and $L$ as

$\displaystyle B_{4}(t,0,L^{-1})$
$\displaystyle=\frac{\partial_{h}^{4}F(t,0,L^{-1})}{(\partial_{h}^{2}F(t,0,L^{-1}))^{2}}+3$
$\displaystyle=\frac{\partial_{h}^{4}\tilde{F}(tL^{y_{t}},0)}{(\partial_{h}^{2}\tilde{F}(tL^{y_{t}},0))^{2}}+3$
$\displaystyle=b_{4}+ctL^{1/\nu}+{\cal O}(t^{2}),$ (8)

for small $t$, where
$b_{4}=\partial_{h}^{4}\tilde{F}(0,0)/(\partial_{h}^{2}\tilde{F}(0,0))^{2}+3$,
$\nu=1/y_{t}$, and $c$ is a constant. Eq. (8) shows that $B_{4}(t,0,L^{-1})$
obtained for various $L$ at $h=0$ has a crossing at $t=0$. The parameter
$b_{4}$ is given only from $\tilde{F}(t,h)$ and thus are specific for each
universality class. For the $Z(2)$ universality class, the value is known to
be [32]

$\displaystyle b_{4}=1.604.$ (9)

Equation (4) means that $F(t,h,L^{-1})$ is the cumulant generating function of
$M$ up to an additive constant. Then, as shown in Eqs. (48) and (50) in
Appendix A, this function is related to the probability distribution function
$p_{M}(M;t,h,L^{-1})$ of $M$ as

$\displaystyle e^{F(t,h^{\prime}-h,L^{-1})}=c_{F}\int
dM\,e^{h^{\prime}M}p_{M}(M;t,h,L^{-1}),$ (10)

where $c_{F}$ is a constant determined from the normalization condition $\int
dM\,p_{M}(M)=1$. Here, let us define another probability distribution
$\tilde{p}_{M}(M;t)$ as

$\displaystyle e^{\tilde{F}(t,h^{\prime})}=c_{\tilde{F}}\int
dMe^{h^{\prime}m}\tilde{p}_{M}(M;t).$ (11)

From Eq. (3), one finds

$\displaystyle
p_{M}(M;t,0,L^{-1})=L^{y_{h}}\tilde{p}_{M}(ML^{-y_{h}};tL^{y_{t}}).$ (12)

When we consider magnetization per unit volume $m=M/V$, the probability
distribution of $m$ is given by

$\displaystyle
p_{m}(m;t,0,L^{-1})=L^{y_{h}-3}\tilde{p}_{M}(mL^{3-y_{h}};tL^{y_{t}}).$ (13)

At the CP, $(t,h)=(0,0)$, one finds from Eq. (13) that

$\displaystyle p_{m}(m;0,0,L^{-1})=L^{y_{h}-3}\tilde{p}_{M}(mL^{3-y_{h}};0).$
(14)

Equation (13) also suggests that, when $\tilde{p}_{M}(M;t)$ has a local
extremum at $M=\tilde{M}(t)$, $p_{m}(m;t,0,L^{-1})$ has corresponding local
extremum at

$\displaystyle m=L^{y_{h}-3}\tilde{M}(tL^{y_{t}}).$ (15)

This implies that the $t$ and $L$ dependences of the maximum of
$p_{m}(m;t,0,L^{-1})$ are described by a single function $\tilde{M}(t)$.

##  III Setup

###  III.1 Lattice action and basic observables

In this study we investigate the four-dimensional system described by the
lattice action of QCD

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

with $S_{g}$ and $S_{q}$ being the gauge and quark actions. For $S_{g}$ we
employ the plaquette action

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

with the gauge coupling parameter $\beta=6/g^{2}$ and the space-time lattice
volume $N_{\rm site}=N_{s}^{3}\times N_{t}$. The plaquette operator $\hat{P}$
is given by

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

where $U_{x,\mu}$ is the link variable in the $\mu$ direction at site $x$,
$x+\hat{\mu}$ is the next site in the $\mu$ direction from $x$, $N_{\rm c}=3$,
and ${\rm tr_{C}}$ is the trace over color indices.

For $S_{q}$, we adopt the Wilson quark action

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

with the Wilson quark kernel

$\displaystyle M_{xy}(\kappa_{f})=\delta_{xy}-\kappa_{f}B_{xy},$ (20)
$\displaystyle
B_{xy}=\sum_{\mu=1}^{4}\left[(1-\gamma_{\mu})\,U_{x,\mu}\,\delta_{y,x+\hat{\mu}}+(1+\gamma_{\mu})\,U_{y,\mu}^{\dagger}\,\delta_{y,x-\hat{\mu}}\right],$
(21)

where $x$, $y$ represent lattice sites. The color and Dirac-spinor indices are
suppressed for simplicity. $\kappa_{f}$ is the hopping parameter for the $f$th
flavor. The bare quark mass $m_{f}$ is related to $\kappa_{f}$ as

$\displaystyle\kappa_{f}=\frac{1}{2am_{f}+8},$ (22)

with the lattice spacing $a$. The matrix $B_{xy}$ has nonzero values only when
lattice sites $x$ and $y$ are located in an adjacent sites. Therefore, this
term represents the “hopping” of a quark between adjacent sites. The heavy
quark limit $m_{f}\to\infty$ corresponds to $\kappa_{f}\to 0$.

In the following, we consider degenerated $N_{\rm f}$ flavors with a common
hopping parameter $\kappa=\kappa_{f}$ corresponding to a common quark mass
$m_{q}$ — generalization to non-degenerate cases is straightforward. In this
case, the expectation value of a gauge operator $\hat{O}(U)$ is calculated as

$\displaystyle\langle\hat{O}(U)\rangle$ $\displaystyle=\frac{1}{Z}\int{\cal
D}U{\cal D}\psi{\cal D}\bar{\psi}\,\hat{O}(U)\,e^{-S_{g}-S_{q}}$
$\displaystyle=\frac{1}{Z}\int{\cal D}U\,\hat{O}(U)[{\rm
det}M(\kappa)]^{N_{\rm f}}e^{-S_{g}}$ $\displaystyle=\frac{1}{Z}\int{\cal
D}U\,\hat{O}(U)\,e^{-S_{g}+N_{\rm f}\ln{\rm det}M(\kappa)},$ (23)

with the partition function $Z=\int{\cal D}U\,e^{-S_{g}+N_{\rm f}\ln{\rm
det}M(\kappa)}$.

In the heavy quark limit $\kappa=0$ ($m_{q}=\infty$), the deconfinement phase
transition at nonzero temperature is characterized by the spontaneous symmetry
breaking of the global $Z(3)$ center symmetry of the $SU(3)$ gauge symmetry.
The most conventional choice for the order parameter of this phase transition
is the Polyakov loop

$\displaystyle\hat{\Omega}=\frac{1}{N_{\rm
c}N_{s}^{3}}\displaystyle\sum_{\vec{x}}{\rm
tr_{C}}\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],$ (24)

where the summation $\sum_{\vec{x}}$ is over the spatial lattice sites on one
time slice. In the heavy quark limit, $\langle\hat{\Omega}\rangle=0$ below the
critical temperature $T_{\rm c}$, while $\langle\hat{\Omega}\rangle$ takes a
nonzero value at $T>T_{\rm c}$. For finite $m_{q}$, the $Z(3)$ symmetry is
explicitly broken by the quark term, and thus $\langle\hat{\Omega}\rangle$
becomes non vanishing for all $T$. Even in this case, when $m_{q}$ is
sufficiently large, $\langle\hat{\Omega}\rangle$ jumps discontinuously at the
first-order transition and thus can be used to detect the first-order
transition line and its CP [27, 29, 30]. In the present study, we focus on the
real part of the Polyakov loop

$\displaystyle\hat{\Omega}_{\rm R}={\rm Re}\,\hat{\Omega},$ (25)

and study its probability distribution function [27, 29]

$\displaystyle p(\Omega_{\rm R})=\langle\delta(\Omega_{\rm
R}-\hat{\Omega}_{\rm R})\rangle.$ (26)

###  III.2 Hopping parameter expansion

Figure 2: Six-step Wilson loops; rectangle (a), chair-type (b) and crown-type
(c). Figure 3: Bent Polyakov loops, $\Omega_{1}$ (left) and $\Omega_{2}$
(middle) for $N_{t}=4$.

To calculate Eq. (23) around the heavy-quark limit, in the present study we
adopt the hopping parameter expansion (HPE) for $\ln{\rm det}M(\kappa)$:

$\displaystyle\ln\left[\frac{\det M(\kappa)}{\det
M(0)}\right]=-\sum_{n=1}^{\infty}\frac{1}{n}{\rm
Tr}\left[B^{n}\right]\kappa^{n},$ (27)

where the matrix $B_{xy}$ is defined by Eq. (21) and ${\rm Tr}$ is the trace
over all indices. In Eq. (27), the contribution at $\kappa=0$ is subtracted
for convenience. Since $B_{xy}$ takes nonzero values only for adjacent lattice
sites $x$ and $y$, $B^{n}$ is graphically represented by trajectories of $n$
links [35]. Because of the trace in Eq. (27), non-vanishing contributions are
given by closed trajectories. The lowest-order contributions of Eq. (27) start
from $n=4$, and all contributions for odd $n$ vanishes when $N_{t}$ is even.
By writing

$\displaystyle-\ln\left[\frac{\det M(\kappa)}{\det M(0)}\right]=$
$\displaystyle S_{\rm LO}+S_{\rm NLO}+{\cal O}(\kappa^{8}),$ (28)

one finds for $N_{t}=4$ that the LO contribution is given by the plaquette and
the Polyakov loop as

$\displaystyle S_{\rm LO}=-2N_{\rm c}\big{(}48N_{\rm
site}\hat{P}+32N_{s}^{3}\hat{\Omega}_{\rm R}\big{)}\kappa^{4}.$ (29)

The NLO term $S_{\rm NLO}$ consists of the six-step Wilson loops and bent
Polyakov loops as

$\displaystyle S_{\rm NLO}=$ $\displaystyle-2N_{\rm c}\big{(}384\,\hat{W}_{\rm
rec}+768\,\hat{W}_{\rm chair}+256\,\hat{W}_{\rm crown}$
$\displaystyle+192\,{\rm Re}\,\hat{\Omega}_{1}+96\,{\rm
Re}\,\hat{\Omega}_{2}\big{)}N_{\rm site}\kappa^{6}.$ (30)

Here, $\hat{W}_{\mathrm{rec}}$, $\hat{W}_{\mathrm{chair}}$, and
$\hat{W}_{\mathrm{crown}}$ represent the six-step Wilson loops of the
rectangular, chair, and crown types, respectively, as illustrated in Fig. 2.
$\hat{\Omega}_{n}$ are the bent Polyakov loops illustrated in Fig. 3, which
run one step in a spatial direction, $n$ steps in the temporal direction and
return to the original line. All the Wilson loops and Polyakov-loop-type loops
are normalized such that
$\hat{W}_{\mathrm{rec}}=\hat{W}_{\mathrm{chair}}=\hat{W}_{\mathrm{crown}}=1$
and $\hat{\Omega}_{n}=1$ in the weak coupling limit, $U_{x,\mu}=1$. Explicit
definitions of these operators as well as the derivation of Eqs. (29) and (30)
are given in Appendix B.

###  III.3 Numerical implementation with HPE

In this study, we generate the gauge configurations with respect to the action
at the LO in the HPE, i.e.

$\displaystyle S_{g+\rm LO}$ $\displaystyle=S_{g}+N_{\rm f}S_{\rm LO}$
$\displaystyle=-2N_{\rm c}N_{\rm site}\big{(}\beta+48N_{\rm
f}\kappa^{4}\big{)}\hat{P}-64N_{\rm c}N_{\rm
f}N_{s}^{3}\kappa^{4}\hat{\Omega}_{\rm R}$ $\displaystyle=-2N_{\rm c}N_{\rm
site}\beta^{*}\hat{P}-\lambda N_{s}^{3}\hat{\Omega}_{\rm R},$ (31)

with

$\displaystyle\beta^{*}=\beta+48N_{\rm f}\kappa^{4},\qquad\lambda=64N_{\rm
c}N_{\rm f}\kappa^{4}.$ (32)

We then perform the measurements at the NLO by incorporating the effect of
$S_{\rm NLO}$ by the multipoint reweighting method [33, 29, 36]. In this
subsection we discuss the numerical implementation of these analyses.

(a) $X_{\mu}(x)$ | (b) $Y_{0}(x)$
---|---
|
Figure 4:  Staple $X_{\mu}(x)$ and the operator $Y_{0}(x)$ in Eq. (36)
corresponding to the link variable $U_{\mu}(x)$ shown by the blue line.

In the Monte Carlo simulations of pure gauge theory, thanks to the locality of
the action $S_{g}$, it is possible to adopt the pseudo-heat-bath (PHB) and
over-relaxation (OR) algorithms for updating gauge configurations. Focusing on
a link variable $U_{\mu}(x)$, the dependence of $S_{g}$ on $U_{\mu}(x)$ is
given by

$\displaystyle\Delta S_{g}(U_{\mu}(x))=-\frac{\beta}{N_{\rm c}}{\rm Re}\ {\rm
tr_{C}}[U_{\mu}(x)X_{\mu}(x)],$ (33)

where the staple

$\displaystyle X_{\mu}(x)=\sum_{\nu\neq\mu}\sum_{s=\pm
1}U_{x+\hat{\mu},s\nu}U_{x+s\hat{\nu},\mu}^{\dagger}U_{x,s\nu}^{\dagger},$
(34)

with $U_{x,-\mu}=U^{\dagger}_{x-\hat{\mu},\mu}$ for $\mu>0$, is graphically
shown in Fig. 4 (a). In the PHB and OR algorithms, the link variable
$U_{\mu}(x)$ is updated according to the probability determined by Eq. (33).
The fact that Eq. (33) is represented only by local variables near
$U_{\mu}(x)$ enables this procedure efficient especially on the memory-
distributed parallel computing.

When the LO term, Eq. (29), is included into the action, the contribution of a
temporal link $U_{0}(x)$ to $S_{g+\rm LO}$ is modified as

$\displaystyle\Delta S_{g+\rm LO}(U_{\mu}(x))$
$\displaystyle=-\frac{\beta^{*}}{N_{\rm c}}{\rm Re}\ {\rm
tr_{C}}\Big{[}U_{0}(x)\Big{(}X_{0}(x)+\frac{\lambda}{\beta^{*}}Y_{0}(x)\Big{)}\Big{]}$
(35)

with

$\displaystyle
Y_{0}(x)=U_{0}(x+\hat{0})U_{0}(x+2\cdot\hat{0})U_{0}(x+3\cdot\hat{0}),$ (36)

which is schematically shown in Fig. 4 (b). For $N_{t}>4$, $Y_{0}(x)$ is given
by the product of $N_{t}-1$ link variables along the temporal direction. The
contribution of a spatial link to $S_{g+\rm LO}$ is unchanged from Eq. (33).

These results on $\Delta S_{g+\rm LO}(U_{\mu}(x))$ suggest that the Monte
Carlo updates of $U_{\mu}(x)$ can be performed by the PHB and OR efficiently
even for $S_{g+\rm LO}$, provided that the temporal direction is not separated
into different parallel nodes and $Y_{0}(x)$ can be calculated efficiently.
Satisfying this condition is not difficult to attain for large-volume
simulations. By taking this advantage, in this study we perform update of
gauge fields at the LO by combining PHB and OR111 After finishing our
numerical analyses, we knew that a similar idea is suggested in Ref. [37]. We
thank F. Karsch for notifying this literature. . Compared with the pure-gauge
simulation, the increase of the numerical cost to deal with $S_{g+\rm LO}$ in
this method is small since the additional multiplications of SU(3) matrices
required for an update are only $N_{t}-2$ times for the temporal links and the
cost to update the spatial links is unchanged.

In the measurement of observables, we include all the contribution of $S_{\rm
NLO}$ using the multipoint reweighting method222 In Ref. [30], the effect of
$S_{\rm NLO}$ is included in part effectively by the effective NLO method. In
this study, we deal with it exactly.. The expectation value of a gauge
observable $\hat{O}(U)$ at the NLO at the parameter set $(\beta,\kappa)$ is
given from the LO simulations at $(\tilde{\beta},\tilde{\kappa})$ as

$\displaystyle\langle\hat{O}(U)\rangle^{\rm NLO}_{\beta,\kappa}$
$\displaystyle=\frac{\int{\cal D}U\,\hat{O}(U)\,e^{-S_{g+\rm
LO}(\beta,\kappa)-S_{\rm NLO}(\beta,\kappa)}}{\int{\cal D}U\,e^{-S_{g+\rm
LO}(\beta,\kappa)-S_{\rm NLO}(\beta,\kappa)}}$ $\displaystyle=\frac{\int{\cal
D}U\,\hat{O}(U)\,e^{-\Delta S_{g+\rm LO}-S_{\rm
NLO}(\beta,\kappa)}e^{-S_{g+\rm LO}(\tilde{\beta},\tilde{\kappa})}}{\int{\cal
D}U\,e^{-\Delta S_{g+\rm LO}-S_{\rm NLO}(\beta,\kappa)}e^{-S_{g+\rm
LO}(\tilde{\beta},\tilde{\kappa})}}$
$\displaystyle=\frac{\langle\hat{O}(U)\,e^{-\Delta S_{g+\rm LO}-S_{\rm
NLO}(\beta,\kappa)}\rangle^{\rm LO}_{\tilde{\beta},\tilde{\kappa}}}{\langle
e^{-\Delta S_{g+\rm LO}-S_{\rm NLO}(\beta,\kappa)}\rangle^{\rm
LO}_{\tilde{\beta},\tilde{\kappa}}},$ (37)

with

$\displaystyle\Delta S_{g+\rm LO}=$ $\displaystyle S_{g+\rm
LO}(\beta,\kappa)-S_{g+\rm LO}(\tilde{\beta},\tilde{\kappa})$ $\displaystyle=$
$\displaystyle-2N_{\rm c}N_{\rm site}\big{(}\beta-\tilde{\beta}+N_{\rm
f}(\kappa^{4}-\tilde{\kappa}^{4})\big{)}\hat{P}$ $\displaystyle+64N_{\rm
c}N_{\rm f}N_{s}^{3}(\kappa^{4}-\tilde{\kappa}^{4})\hat{\Omega}_{\rm R},$ (38)

and $\langle\cdot\rangle^{\rm LO}_{\beta,\kappa}$ is the expectation value
taken with the action $S_{g+\rm LO}(\beta,\kappa)$. We generate gauge
configurations for the LO action for several values of
$(\tilde{\beta},\tilde{\kappa})$, and evaluate the expectation values at the
NLO by the multipoint reweighting method.

###  III.4 Simulation parameters

Table 1: Simulation parameters: lattice size $N_{s}^{3}\times N_{t}$, $\beta^{*}$, and $\lambda$. The value of $\kappa$ corresponding to each $\lambda$ is also listed for the case $N_{\rm f}=2$.  lattice size | $\beta^{*}$ | $\lambda$ | $\kappa^{N_{\rm f}=2}$
---|---|---|---
$48^{3}\times 4$ | 5.6869 | 0.004 | 0.0568
| 5.6861 | 0.005 | 0.0601
| 5.6849 | 0.006 | 0.0629
$40^{3}\times 4$, $36^{3}\times 4$ | 5.6885 | 0.003 | 0.0529
| 5.6869 | 0.004 | 0.0568
| 5.6861 | 0.005 | 0.0601
| 5.6849 | 0.006 | 0.0629
| 5.6837 | 0.007 | 0.0653
$32^{3}\times 4$ | 5.6885 | 0.003 | 0.0529
| 5.6865 | 0.004 | 0.0568
| 5.6861 | 0.005 | 0.0601
| 5.6845 | 0.006 | 0.0629
| 5.6837 | 0.007 | 0.0653
$24^{3}\times 4$ | 5.6870 | 0.0038 | 0.0561
| 5.6820 | 0.0077 | 0.0669
| 5.6780 | 0.0115 | 0.0740

In this study, we perform Monte Carlo simulations with fixed temporal lattice
size $N_{t}=4$, while the spatial extent $N_{s}$ is changed from $24$ to $48$.
This allows us to perform simulations with large aspect ratio $N_{s}/N_{t}=LT$
up to $12$, where $L$ is the lattice size along the spatial direction in
physical units. For each $N_{s}$, the gauge configurations are generated for
$3$ to $5$ sets of $(\beta^{*},\lambda)$ shown in Table. 1, which are chosen
so that $\beta^{*}$ is close to the transition line at the LO.

The gauge configurations are updated with the LO action $S_{g+\rm LO}$ using
the PHB and OR algorithms as discussed in Sec. III.3. Gauge configurations are
updated by five OR steps after each PHB step. We measure observables every two
sets of the PHB+OR updates, i.e. totally ten OR steps and two PHB steps. For
all parameters we have performed $6\times 10^{5}$ measurements in this way. In
the following, we set $N_{\rm f}=2$ to show the numerical results unless
otherwise stated.

Figure 5:  Examples of the Monte Carlo history of $\hat{\Omega}_{\rm R}$.
Horizontal axis represents the Monte Carlo time in the unit of measurements
which are made every two sets of PHB+OR updates.

In Monte Carlo simulations near a first-order transition, because the
transition between the coexisting two phases becomes rare when the spatial
volume of the system is large, observables averaged over the two phases tend
to have quite long autocorrelations. In Fig. 5 we show the Monte Carlo time
history of $\hat{\Omega}_{\rm R}$ for $N_{s}/N_{t}=LT=10$ and $12$ at the
smallest $\lambda$, at which the autocorrelation is the longest. The
horizontal axis represents the Monte Carlo time in the unit of measurements.
The figure shows that flip-flops between the two phases occur several times in
$10,000$ measurements. The autocorrelation lengths estimated from the
correlation function are about $1,900$ and $800$ measurements for these time
histories, respectively. Throughout this study, we estimate the statistical
errors of observables by the jackknife method unless otherwise stated,
adopting the binsize of $10,000$ measurements which is sufficiently larger
than the estimated autocorrelation lengths. We checked that the statistical
errors thus estimated are roughly stable within a variation of the binsize
from $5,000$ to $30,000$ measurements.

Figure 6:  Probability distribution $p(P,\Omega_{\rm R})$ of the gauge
configurations on the $40^{3}\times 4$ lattice for the simulation points at
$\lambda=0.003$ and $0.007$. Solid lines represent the contours of
$p(P,\Omega_{\rm R})$ obtained with the LO action. The label on each contour
shows the probability inside the contour. Dotted lines represents the contour
lines of $p(P,\Omega_{\rm R})$ with the NLO action at the same parameters
obtained by reweighting the LO data.

###  III.5 Overlapping problem

Our main objective to use the LO action for configuration generation is to
avoid the overlapping problem in the reweighting. In Fig. 6 we show the
contour plot for the probability distribution function of $\hat{P}$ and
$\hat{\Omega}_{\rm R}$,

$\displaystyle p(P,\Omega_{\rm
R})=\left\langle\delta(P-\hat{P})\delta(\Omega_{\rm R}-\hat{\Omega}_{\rm
R})\right\rangle,$ (39)

obtained at $\lambda=0.003$ (red) and $0.007$ (blue) on $40^{3}\times 4$
lattices. The solid lines are the contours for the distribution measured on
the LO configurations generated with $S_{g+\rm LO}$. Each contour curve is
drawn such that the probability inside the contour is $0.9$, $0.7$, $\cdots$,
and $0.1$. To smoothen the plot, we approximate the delta function in Eq. (39)
by the Gauss function,
$\delta(x)\simeq\exp[-(x/\Delta)^{2}]/(\Delta\sqrt{\pi})$, with the width
along $\Omega_{\rm R}$ and $P$ directions $\Delta_{\Omega_{\rm R}}=0.002$ and
$\Delta_{P}=0.0001$, respectively.

In Fig. 6, we also show by the dotted lines the contours of $p(P,\Omega_{\rm
R})$ at the NLO, which is obtained by reweighting the LO data at the same
$(\beta^{*},\lambda)$. The meaning of the contours is the same as the solid
lines. We find that the deviation of the NLO distribution from the original
one at the LO is not significant, suggesting that the effects of the NLO
contribution are not large. The large overlap of the LO and NLO distributions
ensures that, for observables constructed from $\hat{P}$ and
$\hat{\Omega}_{\rm R}$, the NLO results obtained by reweighting the LO data at
the same $\beta^{*}$ and $\lambda$ are statistically reliable. In the analysis
of the CP, the overlapping of the distributions is even more improved after
adjusting the parameters to the transition line.

On the other hand, from this figure, we find that the overlapping of the
distributions at $\lambda=0.003$ and $0.007$ is quite poor — the regions with
probability larger than 0.5 are not overlapping at all with each other. This
means that, if we were to calculate observables at $\lambda=0.007$ by
reweighting data obtained at $\lambda=0.003$, or vice versa, the statistical
quality of the results would be quite low. In Refs. [27, 29, 30], the CP in
heavy quark region was investigated on lattices with $N_{s}/N_{t}=4$–$6$, by
reweighting from pure gauge configurations, i.e. those obtained at
$\lambda=0$. Because the overlapping problem becomes quickly severe as the
system volume becomes large, the same strategy is not applicable to the
present study in which much larger system volumes up to $N_{s}/N_{t}=LT=12$
are simulated. Figure 6 shows that, to incorporate the NLO effects by
reweighting, the use of the LO action for configuration generation is
sufficiently effective in suppressing the overlapping problem. The smallness
of the NLO effects in Fig. 6 further suggests that the effects of dynamical
quarks are dominated by the LO term for these values of $\lambda$.

##  IV Binder cumulant analysis

###  IV.1 Transition line

Figure 7:  Transition line $\beta^{*}_{\rm tr}(\lambda)$ as function of the
aspect ratio $LT=N_{s}/N_{t}$, determined through three definitions: (1)
minimum of $B_{4}^{\Omega}$, (2) zero of $\langle\Omega_{\rm
R}^{3}\rangle_{\rm c}$, and (3) maximum of $\langle\Omega_{\rm
R}^{2}\rangle_{\rm c}$. The results for $\lambda=0.003$, 0.005, 0.007 are
shown from left to right.

We first determine the location of the transition line that corresponds to
$h=0$ in terms of the Ising parameters; see Fig. 1. In the coupling parameter
space $(\beta^{*},\lambda)$, we denote the transition line as $\beta^{*}_{\rm
tr}(\lambda)$. In this study, we determine $\beta^{*}_{\rm tr}$ at each
$\lambda$ adopting the following three conventional choices:

  * •

Maximum of $\langle\Omega_{\rm R}^{2}\rangle_{\rm c}$

  * •

Zero point of $\langle\Omega_{\rm R}^{3}\rangle_{\rm c}$

  * •

Minimum of $B_{4}^{\Omega}=\langle\Omega_{\rm R}^{4}\rangle_{\rm
c}/\langle\Omega_{\rm R}^{2}\rangle_{\rm c}^{2}+3$

In Fig. 7, we show the $LT$ dependence of $\beta^{*}_{\rm tr}$ determined by
these definitions for several values of $\lambda$. The figure shows that the
maximum of $\langle\Omega_{\rm R}^{2}\rangle_{\rm c}$ has a visible $LT$
dependence. On the other hand, the zero point of $\langle\Omega_{\rm
R}^{3}\rangle_{\rm c}$ and the minimum of $B_{4}^{\Omega}$ do not have
statistically-significant $LT$ dependence for $LT\geq 8$. This result shows
that the zero point of $\langle\Omega_{\rm R}^{3}\rangle_{\rm c}$ and minimum
of $B_{4}^{\Omega}$ are sufficiently close to the $\beta^{*}_{\rm tr}$ in the
$L\to\infty$ limit in this range of $LT$. In the following, we employ the
minimum of $B_{4}^{\Omega}$ for the definition of the transition line
$\beta^{*}_{\rm tr}$ for each $N_{t}$. In Fig. 8, we show the transition line
on the $(\beta^{*},\lambda)$ and $(\beta,\lambda)$ planes obtained at $LT=10$.

Figure 8:  Transition line $\beta^{*}_{\rm tr}$ on the $(\beta,\lambda)$ and
$(\beta^{*},\lambda)$ plane obtained at $LT=10$. The points with the error bar
are the results obtained on each measurement. The shaded areas show the result
obtained by the multipoint reweighting method with the error band representing
the statistical error. The dotted vertical line shows the value of
$\lambda_{\rm c}$ determined by the analysis in Sec. IV.2.

Figure 9:  Distribution of $\Omega_{\rm R}$ at $\lambda=0.003$, $0.005$ and
$0.007$.

In Fig. 9, we show the distribution function of $\Omega_{\rm R}$, Eq. (26), on
the transition line for several values of $\lambda$, where the delta function
in Eq. (26) is smeared by the Gaussian approximation as before with the width
$\Delta_{\Omega_{\rm R}}=0.002$. The shaded bands represent the statistical
errors. At $\lambda=0.003$, we see a clear two-peak structure in
$p(\Omega_{\rm R})$ and find that the peaks become sharper as $LT$ becomes
larger. This behavior suggests the first-order phase transition at
$\lambda=0.003$. At $\lambda=0.007$, on the other hand, while two peaks are
observed for $LT\leq 9$, the two peaks cease to exist as $LT$ becomes large.
This suggests the crossover transition in the $L\to\infty$ limit at this
$\lambda$.

###  IV.2 Binder cumulant

Figure 10:  Binder cumulant $B_{4}^{\Omega}$ as a function of $\lambda$
obtained at five $LT=N_{s}/N_{t}$. The statistical errors are shown by the
shaded area. The lower panel is an enlargement of the upper panel around the
crossing point, where the dotted rectangle in the upper panel represents the
region of the lower panel. The points in the lower panel with error bars show
the results of the four parameter fit with Eq. (41). See text for details.

Next, let us determine the position of the CP on the $(\beta,\kappa)$ plane.
As discussed in Sec. II, it is convenient to employ the Binder cumulant
$B_{4}$ of $\Omega_{\rm R}$

$\displaystyle B_{4}^{\Omega}=\frac{\langle\Omega_{\rm R}^{4}\rangle_{\rm
c}}{\langle\Omega_{\rm R}^{2}\rangle_{\rm c}^{2}}+3.$ (40)

This quantity approaches the known values given in Eq. (7) in the $L\to\infty$
limit depending on the order of the transition. Furthermore, provided that
$\Omega_{\rm R}$ corresponds to $m=M/V$ of the Ising model, $B_{4}^{\Omega}$
should obey Eq. (8) near the CP.

In the upper panel of Fig. 10, we show $B_{4}^{\Omega}$ along the transition
line as a function of $\lambda$ for five values of $LT$. $\lambda$ is varied
continuously by the multipoint reweighting method. The lower panel is an
enlargement of the upper panel around the crossing point. The figure shows
that $B_{4}^{\Omega}$ has a crossing at $\lambda=\lambda_{\rm c}\simeq 0.005$
and is an increasing (decreasing) function of $LT$ for $\lambda>\lambda_{\rm
c}$ ($\lambda<\lambda_{\rm c}$). The existence of the CP at $\lambda\simeq
0.005$ is suggested from this result.

Table 2:  Results of the four parameter fits of $B_{4}^{\Omega}$ with Eq. (41). Three largest volumes with $LT=N_{s}/N_{t}=12,~{}10,~{}9$ are used for the fits. The left column shows the values of $\lambda$ used for the fit.  fit $\lambda~{}(\times 10^{4})$ | $b_{4}$ | $\lambda_{\rm c}$ | $\nu$ | $\chi^{2}/{\rm dof}$
---|---|---|---|---
49, 52 | 1.631(24) | 0.00503(14) | 0.614(49) | 0.45
48, 53 | 1.630(24) | 0.00503(14) | 0.614(48) | 0.46
45, 55 | 1.629(24) | 0.00502(14) | 0.622(48) | 0.46
45, 50 | 1.630(24) | 0.00503(15) | 0.634(47) | 0.48
50, 55 | 1.631(24) | 0.00504(14) | 0.610(51) | 0.38
45, 50, 55 | 1.620(23) | 0.00494(14) | 0.626(41) | 7.9

To determine $\lambda_{\rm c}$ and the critical exponent $\nu$ quantitatively,
we fit the numerical results of $B_{4}^{\Omega}$ by a fitting function
motivated by Eq. (8):

$\displaystyle B_{4}^{\Omega}(\lambda,LT)=b_{4}+c(\lambda-\lambda_{\rm
c})(LT)^{1/\nu},$ (41)

where $b_{4}$, $\lambda_{\rm c}$, $\nu$, $c$ are the fit parameters. In this
study, we can vary $\lambda$ continuously by the multipoint reweighting
method. However, because data at different $\lambda$ on the same volume are
correlated, it is not meaningful to use too many $\lambda$ values. Using the
data at three largest volumes, $LT=12,~{}10,~{}9$, two or three $\lambda$
values (6 or 9 data points, respectively) should be sufficient for the four
parameter fit of Eq. (41). We thus repeat the fit for several choices of
$\lambda$ values, taking the covariance between data at different $\lambda$
into account in the calculation of $\chi^{2}$.

In Table 2, we summarize the results of the fit using the data at three
largest volumes, $LT=12,~{}10,~{}9$, and at $\lambda$ values listed in the
left column of the table. The statistical error in the table is estimated by
the standard chi-square analysis. The table shows that the value of
$\chi^{2}/{\rm dof}$ are smaller than unity in the fits with two $\lambda$
values, but $\chi^{2}/{\rm dof}$ is unacceptably large with three $\lambda$
values, while all the results for the fitting parameters are well consistent
within errors. We choose the fit result for $\lambda=(0.0048,~{}0.0053)$
depicted by bold characters in Table 2 as the central value and include the
uncertainty in the fits with two $\lambda$ values as the systematic error.

Table 3:  Results of the four parameter fit of $B_{4}^{\Omega}$ using the data points at three and four largest and all volumes. The first parentheses are for statistical errors estimated by the jackknife method, and the second parentheses are for systematic errors due to the choice of $\lambda$ values for the fit. See Sec. IV.2 for details.  $LT=N_{s}/N_{t}$ | $b_{4}$ | $\lambda_{\rm c}$ | $\nu$ | $\chi^{2}/{\rm dof}$
---|---|---|---|---
12, 10, 9 | 1.630(24)(2) | 0.00503(14)(2) | 0.614(48)(3) | 0.46
12, 10, 9, 8 | 1.643(15)(2) | 0.00510(10)(2) | 0.614(29)(3) | 0.37
12, 10, 9, 8, 6 | 1.645(11)(2) | 0.00511(8)(2) | 0.593(18)(3) | 0.67

We repeat the analysis also with other sets of system volumes. The results of
the fits with four and five largest volumes, together with the fit with three
largest volumes, are summarized in Table 3. For the fit with four and five
largest volumes, $LT=12$–$8$ and $12$–$6$, we now have 8 and 10 data points
for the four parameter fit with two $\lambda$ values. For these fits, we
choose $\lambda=(0.0048,0.0053)$ as the central value again.

The results of $b_{4}$ and $\lambda_{\rm c}$ obtained by these fits are shown
in the lower panel of Fig. 10. The thick symbols show the central value, and
thin symbols with light colors show the results obtained by the variation of
$\lambda$ in Table 2. The results with the three largest volumes are shown by
black triangles, while those with four and five largest volumes are shown by
blue squares and green pentagons, respectively. In the figure, $b_{4}$
expected from the $Z(2)$ universality, Eq. (9), is shown by the dashed
horizontal line.

From Table 3 and Fig. 10, we find that, when we adopt the fit with the three
largest volumes, $LT\geq 9$, the fit result of $b_{4}$ is consistent with the
$Z(2)$ value within about $1\sigma$. On the other hand, when we include
smaller volumes, $LT\geq 8$ or $LT\geq 6$, $b_{4}$ from the fits show
statistically significant deviation from the $Z(2)$ value. In Table 3, we also
summarize the results for the critical exponent $\nu$. From the $Z(2)$
universality class, we expect $\nu=1/y_{t}\approx 0.630$. We find that the
result of the fit with $LT\geq 9$ is consistent with the $Z(2)$ value within
the error, while the result of the fit with $LT\geq 6$ has a significant
deviation from the $Z(2)$ value, though the values of $\chi^{2}/{\rm dof}$ are
all smaller than unity.

We thus conclude that the FSS in the $Z(2)$ universality class is confirmed
when the system volume is large enough, $LT\geq 9$ — lattices with $LT\leq 8$
are not large enough to apply a FSS analysis for $B_{4}^{\Omega}$. The value
of $\lambda_{\rm c}$ thus determined is also shown in Fig. 8.

In Appendix C, we perform the analysis of $B_{4}^{\Omega}$ at the LO of the
HPE and compare the results with those at the NLO discussed in this Section.
We find that the LO result for $\lambda_{\rm c}$ is about 2.6% larger than the
NLO value. This small difference suggests that the truncation error of the HPE
is well under control at the NLO around $\lambda_{\rm c}$.

###  IV.3 Mixing with energy-like observable

So far, we have performed the analyses of $B_{4}^{\Omega}$ assuming that
$\Omega_{\rm R}$ corresponds to the magnetization $m=M/V$ of the Ising model.
Although our numerical results thus far are in good agreement with this
assumption, a possible mixing with the energy-like observable [14, 22] in
$\Omega_{\rm R}$ is not excluded in general. In this case, the behavior of
$B_{4}^{\Omega}$ near the CP is modified from Eq. (8) as [22]

$\displaystyle B_{4}^{\Omega}(\lambda,LT)=$
$\displaystyle\big{(}b_{4}+c(\lambda-\lambda_{\rm
c})(LT)^{1/\nu}\big{)}\big{(}1+d(LT)^{y_{t}-y_{h}}\big{)}.$ (42)

To investigate the effect of this possible mixing, we try fits of
$B_{4}^{\Omega}$ based on Eq. (42). We use the values of $B_{4}^{\Omega}$ at
three $\lambda$ for the fits to increase the number of data points. We find
that the six parameter fits with Eq. (42) with the fitting parameters $b_{4}$,
$\lambda_{\rm c}$, $\nu$, $c$, $d$, $y_{t}-y_{h}$ are quite unstable,
suggesting that $\chi^{2}$ has many local minima. The model space of Eq. (42)
would be too large against the data. As a next trial, we perform five
parameter fits with Eq. (42) by fixing $y_{t}-y_{h}=-0.894$. In this case, we
find that $\chi^{2}$ still has many local minima, and $\chi^{2}/{\rm dof}$
becomes larger compared with the four parameter fit. It is also found that the
value of $d$ is consistent with zero within the error for all trials with the
variation of $\lambda$ values. This suggests that the mixing of the energy-
like observable in $\Omega_{\rm R}$ is negligible around the CP in the heavy-
quark region.

Table 4:  Location of the critical point $(\beta_{\rm c},\kappa_{\rm c})$ for various $N_{\rm f}$. For $\lambda_{\rm c}$, the first parentheses are for statistical errors and the second parentheses are for systematic errors from the fit as discussed in Sec. IV.2. The errors for $\beta_{\rm c}$ and $\kappa_{\rm c}$ include the systematic errors.  $N_{\rm f}$ | $\beta_{\rm c}$ | $\kappa_{\rm c}$ | $\lambda_{\rm c}$
---|---|---|---
$1$ | 5.68446(22) | 0.0714(5) | 0.00498(14)(2)
$2$ | 5.68453(22) | 0.0602(4) | 0.00503(14)(2)
$3$ | 5.68456(21) | 0.0544(4) | 0.00505(14)(2)

###  IV.4 $N_{\rm f}$ dependence

In Table 4, we summarize our final results for the location of the CP,
$(\beta_{\rm c},\kappa_{\rm c})$. In the table, we also show the results for
$N_{\rm f}=1$ and $3$. We note that the $N_{\rm f}$ dependence of the HPE is
trivial at the LO in the sense that $N_{\rm f}$ enters the action Eq. (31) at
this order only through the combination $\lambda=64N_{\rm c}N_{\rm
f}\kappa^{4}$ after the replacement $\beta\to\beta^{*}$. Therefore,
$\lambda_{\rm c}$ does not depend on $N_{\rm f}$. At the LO, this allows us to
obtain the value of $\kappa_{\rm c}$ for various $N_{\rm f}$ from the value of
$\kappa_{\rm c}$ at $N_{\rm f}=2$ [29]. Because such a simple scaling is no
longer applicable at the NLO, we made individual numerical analyses at $N_{\rm
f}=1$ and $3$. From Table 4, we find that the results of $\lambda_{\rm c}$ are
almost insensitive to $N_{\rm f}$. This means that the NLO effects on
$\lambda_{\rm c}$ are small.

##  V Distribution function of $\Omega_{\rm R}$

In this section, we study the scaling behavior of the distribution function
$p(\Omega_{\rm R})$ to further investigate the consistency with the $Z(2)$
universality class around the CP.

###  V.1 Scaling of distribution function

Figure 11:  Effective potential $V(\Omega_{\rm R})=-\ln p(\Omega_{\rm R})$.
Bottom panel is an enlargement of the region enclosed by a dotted rectangle in
the top panel.

Let us first focus on the $LT$ dependence of $p(\Omega_{\rm R})$ at the CP. In
the following, instead of $p(\Omega_{\rm R})$ itself, we study the effective
potential defined from $p(\Omega_{\rm R})$:

$\displaystyle V(\Omega_{\rm R};\lambda,LT)=-\ln p(\Omega_{\rm
R})_{\lambda,LT},$ (43)

as this quantity is more convenient in comparing the results at different $LT$
[27, 29]. From Eq. (14), the $LT$ dependence of $V(\Omega_{\rm R},\lambda,LT)$
at the CP will be described by a single function $\tilde{V}(x)$ as

$\displaystyle V(\Omega_{\rm R};\lambda_{\rm
c},LT)=\tilde{V}\Big{(}\big{(}\Omega_{\rm R}-\langle\Omega_{\rm
R}\rangle\big{)}(LT)^{3-y_{h}}\Big{)},$ (44)

up to an additive constant, where $\langle\Omega_{\rm R}\rangle$ is subtracted
from $\Omega_{\rm R}$ to adjust the center of the distribution.

To see if the scaling behavior of Eq. (44) is satisfied, we show in Fig. 11
the effective potential $V(\Omega_{\rm R};\lambda_{\rm c},LT)$ at the CP
obtained at five values of $LT$, as a function of $\big{(}\Omega_{\rm
R}-\langle\Omega_{\rm R}\rangle\big{)}(LT)^{3-y_{h}}$, where we set
$3-y_{h}=0.504$. For the figure, we adjust the arbitrary constant term of
$V(\Omega_{\rm R};\lambda_{\rm c},LT)$ such that $V(\Omega^{(1)};\lambda_{\rm
c},LT)+V(\Omega^{(2)};\lambda_{\rm c},LT)=0$, where $\Omega^{(1)}$ and
$\Omega^{(2)}~{}(>\Omega^{(1)})$ are the values of $\Omega_{\rm R}$ at the two
local minima of $V(\Omega_{\rm R};\lambda_{\rm c},LT)$. No further adjustments
are made in the figure. The error bands do not include the uncertainty of the
additive constant. The lower panel is an enlargement of the region indicated
by the dotted rectangle in the upper panel. From Fig. 11, we find that the
numerical results for $LT=8$–12 agree almost completely within the errors with
the scaling relation Eq. (44). This result nicely supports the FSS in the
$Z(2)$ universality class at the CP. From the upper panel of Fig. 11, we note
that the effective potential for $LT=6$ shows a clear deviation from the
results for larger volumes at $\Omega_{\rm R}\ll\Omega^{(1)}$ and $\Omega_{\rm
R}\gg\Omega^{(2)}$, while it agrees well with them in the range
$\Omega^{(1)}\lesssim\Omega_{\rm R}\lesssim\Omega^{(2)}$. This suggests that
the deviation from the $Z(2)$ FSS by lattices with small $LT$, discussed in
Sec. IV.2, is due to that in the tails of the distribution $p(\Omega_{\rm R})$
for small $LT$.

###  V.2 Gap between the two minima

Using Eq. (13), the argument of Sec. V.1 on the effective potential can be
extended away from the CP along the transition line. In this subsection, we
study the gap between the two local minima of $V(\Omega_{\rm R};\lambda_{\rm
c},LT)$,

$\displaystyle\Delta\Omega=\Omega^{(2)}-\Omega^{(1)}.$ (45)

According to Eq. (15), this quantity should behave around the CP as

$\displaystyle\Delta\Omega(\lambda,LT)=(LT)^{y_{h}-3}\Delta\tilde{\Omega}\big{(}(\lambda-\lambda_{\rm
c})(LT)^{1/\nu}\big{)}$ (46)

provided that $p(\Omega_{\rm R})$ obeys the FSS.

Figure 12:  Positions of peaks of the distribution function $p(\Omega_{\rm
R})$ measured on the transition line.

In Fig. 12, we show the $\lambda$ dependence of $\Omega^{(1)}$ and
$\Omega^{(2)}$. As seen from Fig. 9, a clear two peak structure of
$p(\Omega_{\rm R})$ disappears when $\lambda$ exceeds some value depending on
$LT$. Even before the disappearance of the two peaks, identification of local
maxima of $p(\Omega_{\rm R})$ becomes unstable by statistical fluctuations. In
Fig. 12, we thus truncate the plots for $\Omega^{(1)}$ and $\Omega^{(2)}$ at
finite $\lambda$. The shaded areas in the figure represent statistical errors
estimated by the jackknife method, for which we repeat the analysis of
$\Omega^{(1,2)}$ for $p(\Omega_{\rm R})$ obtained in each jackknife sample
with the smearing width of $\Delta_{\Omega_{\rm R}}=0.002$. 333 We see that
the errors in Fig. 12 become occasionally large. We find that this is due to
statistical oscillations in the shape of $p(\Omega_{\rm R})$ around the peak:
Though the oscillations are within the statistical errors, the peak position
in each jackknife bin can jump discontinuously when oscillation appears just
at the peak position as we vary $\lambda$. This makes the resulting jackknife
error large there. From this observation, we think that these large errors are
overestimated.

Figure 13:  Scaling of the gap $\Delta\Omega$ around $T_{\rm c}$.

From Fig. 12 we extract $\Delta\Omega$ as a function of $\lambda$. In Fig. 13,
we show $\Delta\Omega$ for five different volumes. To see the FSS, the
vertical and horizontal axes are adjusted according to Eq. (46), where the
$Z(2)$ values $3-y_{h}=0.519$ and $\nu=0.630$, and $\lambda_{\rm c}=0.00503$
determined in the previous section are used. The figure shows that, for a wide
range of $\lambda-\lambda_{\rm c}$ and $LT$, the results of $\Delta\Omega$
obtained on different volumes are on top of each other within the errors. This
supports the FSS of $p(\Omega_{\rm R})$ around the peak positions over a wide
range of $LT$ and $\lambda$.

It is interesting to note that the scaling behavior of $\Delta\Omega$ is
observed even at $LT=6$, although the FSS of $B_{4}^{\Omega}$ is violated
already at $LT=8$. As discussed in the previous subsection, we may understand
this when the violation of the FSS for $B_{4}^{\Omega}$ is due to the
violation in the tails of the distribution function $p(\Omega_{\rm R})$. As
seen in Fig. 11, $V(\Omega_{\rm R})$ at various volumes agrees well for
$\Omega^{(1)}\lesssim\Omega\lesssim\Omega^{(2)}$ even for small values of
$LT$. As the higher-order cumulants are sensitive to the whole structure of
the distribution, $B_{4}^{\Omega}$ will be more sensitive to the violation of
FSS at the tails of the distribution. On the other hand, $\Delta\Omega$ is
insensitive to them by definition. We also note that the statistical error of
$\Delta\Omega$ is naturally small because it is defined by the peaks of the
distribution. Therefore, $\Delta\Omega$ is useful in seeing the FSS around the
CP.

Figure 14:  Gap $\Delta\Omega$ around $T_{\rm c}$ as a function of
$(LT)^{3-y_{h}}$. The $Z(2)$ value $3-y_{h}=0.519$ is assumed.

From Eq. (46), $\Delta\Omega$ should behave linearly as a function of
$(LT)^{3-y_{h}}$ at the CP $\lambda=\lambda_{\rm c}$. In Fig. 14, we show
$\Delta\Omega$ on the transition line at various values of $\lambda$, as a
function of $(LT)^{3-y_{h}}$. In the same figure, the dashed lines show linear
functions $\Delta\Omega=k(LT)^{3-y_{h}}$ for various values of $k$. Figure 14
suggests that the linear behavior is realized at $\lambda\simeq 0.005$, which
is consistent with our estimation $\lambda_{\rm c}=0.00503(14)(2)$ from the
analysis of $B_{4}^{\Omega}$.

###  V.3 Discussions

Let us comment on the relation of the present results with those given in
Refs. [27, 29, 30]. In these studies, the CP is defined as the point at which
the two peak structure of $p(\Omega_{\rm R})$ disappears. On finite lattice
with finite $LT$, this leads to $\lambda$ which is larger than our value of
$\lambda_{\rm c}$ in the $L\to\infty$ limit. In fact, in Refs. [27], the
location of the CP is estimated as $\kappa_{\rm c}=0.0658(3)(^{+4}_{-11})$ for
$N_{\rm f}=2$ and $N_{t}=4$, which is about 10% larger than that given in
Table 4. Values of $\lambda_{\rm c}$ ($\kappa_{\rm c}$) which are smaller than
Ref. [30] for each $N_{t}$ are also reported by a recent study of the CP in
the heavy quark region on fine lattices ($N_{t}=6$–$10$) using the Binder
cumulant method [31]. Though the difference may be removed in the $L\to\infty$
limit, a careful extrapolation will be required. Because the FSS is clearly
identified in this study, we think that the extrapolation to the large volume
limit is stably performed with the present analysis.

We also note that the latent heat at the deconfinement transition in the
$SU(3)$ Yang-Mills theory ($\kappa=0$) has been measured in Ref. [38]
recently. It was found that the latent heat becomes larger with increasing the
spatial volume. This may be attributed to a remnant of the FSS around the
$Z(2)$ CP like $\Delta\Omega$ studied in the present study.

##  VI Conclusions

In this paper, we studied the distribution function of the Polyakov loop and
its cumulants around the CP in the heavy quark region of QCD. Large volume
simulations up to the aspect ratio $N_{s}/N_{t}=LT=12$ have been carried out
to see the finite-size scaling, while the lattice spacing is fixed to
$N_{t}=4$. We have performed the measurement of observables using the hopping
parameter expansion for the quark determinant; the measurement has been
performed at the next-to-leading order of the hopping parameter expansion by
the multipoint reweighting method evaluated on the gauge configurations
generated for the action at the leading order. We found that this analysis is
quite effective in reducing statistical errors by avoiding the overlapping
problem of the reweighting method, while the numerical cost hardly changes
from the pure Yang-Mills simulations. The convergence of the hopping parameter
expansion at the next-to-leading order at the critical point is also verified
by the comparison with the leading order result.

Using the data on $p(\Omega_{\rm R})$ thus obtained, we have performed the
Binder cumulant analysis for determining the location of the critical point
and evaluating the critical exponent. We found that the critical exponent
$\nu$ and the value of the Polyakov-loop Binder cumulant $B_{4}^{\Omega}$ at
the critical point is consistent with the $Z(2)$ universality class when
$LT\geq 9$ data are used for the analysis. On the other hand, statistically-
significant deviation from the $Z(2)$ scaling is observed when the data at
$LT=8$ is included, which suggests that this spatial volume is not large
enough to apply the finite-size scaling.

The scaling behavior near the critical point is further studied using the
structure of the distribution function of the real part of the Polyakov loop,
$p(\Omega_{\rm R})$. We found that the structure of $p(\Omega_{\rm R})$ for
various $LT$ obeys the finite-size scaling well especially near the peaks of
$p(\Omega_{\rm R})$. We have also proposed the use of the gap $\Delta\Omega$
between the peaks of $p(\Omega_{\rm R})$ for the finite-size scaling analysis.
We showed that the $\lambda$ and $LT$ dependence of $\Delta\Omega$ is in good
agreement with the $Z(2)$ scaling over a wide range of $\lambda$ and $LT$. On
the other hand, the deviation of $p(\Omega_{\rm R})$ around the tails of the
distribution is observed on small lattices, which would give rise to the
violation of the finite-size scaling of $B_{4}^{\Omega}$ in small volumes.

Acknowledgments

The authors thank Frithjof Karsch, Makoto Kikuchi, Yoshifumi Nakamura, and the
members of WHOT-QCD Collaboration for useful discussions. This work was
supported by in part JSPS KAKENHI (Grant Nos. JP17K05442, JP18K03607,
JP19H05598, JP19K03819, JP19H05146, JP21K03550), the Uchida Energy Science
Promotion Foundation, the HPCI System Research project (Project ID: hp170208,
hp190036, hp200089, hp210039), and Joint Usage/Research Center for
Interdisciplinary Large-scale Information Infrastructures in Japan (JHPCN)
(Project ID: jh190003, jh190035, jh190063, jh200049). This research used
computational resources of OCTPUS of the large-scale computation program at
the Cybermedia Center, Osaka University, and ITO of the JHPCN Start-Up
Projects at the Research Institute for Information Technology, Kyushu
University.

##  Appendix A Cumulants

Let us consider a probability distribution function $p(x)$ of a stochastic
variable $x$. Since $p(x)$ represents probability, it satisfies the
normalization condition $\int\\!dx\,p(x)=1$.

The $m$th-order moment $\langle x^{m}\rangle$ of $p(x)$ is defined by

$\displaystyle\langle x^{m}\rangle=\int\\!dx\,x^{m}p(x).$ (47)

Using the moment generating function

$\displaystyle G(\theta)=\int\\!dx\,e^{x\theta}p(x)=\langle
e^{x\theta}\rangle,$ (48)

the moments are also given by

$\displaystyle\langle
x^{m}\rangle=\partial_{\theta}^{m}G(\theta)|_{\theta=0},$ (49)

with $\partial_{\theta}=\partial/\partial\theta$.

The cumulants are defined from the cumulant generating function

$\displaystyle K(\theta)=\ln G(\theta)=\ln\langle e^{x\theta}\rangle$ (50)

as

$\displaystyle\langle x^{m}\rangle_{\rm
c}=\partial_{\theta}^{m}K(\theta)|_{\theta=0}.$ (51)

From Eqs. (49) and (51), one easily finds that the cumulants are related to
the moments as

$\displaystyle\langle x\rangle_{\rm c}$ $\displaystyle=\langle x\rangle,$ (52)
$\displaystyle\langle x^{2}\rangle_{\rm c}$ $\displaystyle=\langle
x^{2}\rangle-\langle x\rangle^{2}=\langle\delta x^{2}\rangle,$ (53)
$\displaystyle\langle x^{3}\rangle_{\rm c}$ $\displaystyle=\langle\delta
x^{3}\rangle,$ (54) $\displaystyle\langle x^{4}\rangle_{\rm c}$
$\displaystyle=\langle\delta x^{4}\rangle-3\langle\delta x^{2}\rangle^{2},$
(55)

and etc. with $\delta x=x-\langle x\rangle$. The cumulants are useful in
representing properties of $p(x)$ than the moments for various purposes. In
particular, in statistical mechanics the cumulants of an extensive variable
are extensive variables; see, for example, Ref. [5].

##  Appendix B Hopping parameter expansion

In this appendix we derive Eqs. (29) and (30). Throughout this Appendix we
assume general value for $N_{t}$.

Figure 15:  “Appendix” structure of trajectories.

As in Eq. (27), the HPE of $\ln[\det M(\kappa)]$ is given by ${\rm
Tr}[B^{n}]$. Since the matrix $B$ has nonzero contributions only between
neighboring lattice sites, ${\rm Tr}[B^{n}]$ are graphically represented by
the closed trajectories with $n$ links [35]. However, trajectories including
“appendices” shown in Fig. 15 do not contribute to the HPE because for such a
path the product of the matrix in Dirac space vanishes at the tip of the
appendix as $(1-\gamma^{\mu})(1+\gamma^{\mu})=0$ [35]. With this exception,
all possible closed trajectories composed of $n$ links contribute to the HPE
at the order $\kappa^{n}$. In the following, we calculate their contributions
by classifying the trajectories by the shape.

Let us start from the plaquette, i.e. $1\times 1$ rectangle, which gives a
lowest-order contribution to Eq. (27) at the order $\kappa^{4}$. The plaquette
operator $\hat{P}$ in Eq. (18) is defined in such a way that
$\langle\hat{P}\rangle=1$ in the weak coupling limit with $U_{\mu}(x)=1$. To
satisfy this condition Eq. (18) has a coefficient $1/(N_{\rm c}M_{\rm
plaq})=1/18$, where $M_{\rm plaq}$ is the number of different plaquettes per
lattice site; $M_{\rm plaq}={}_{4}{\rm C}_{2}=6$, which is the number of
combinations of axes $(\mu,\nu)$ at which the plaquettes are located.

The contribution from all plaquettes to Eq. (27) is calculated to be

$\displaystyle-2N_{\rm c}M_{\rm plaq}D_{\rm plaq}N_{\rm
site}\kappa^{4}\hat{P}.$ (56)

Here, $D_{\rm plaq}$ is the coefficient from the trace in the Dirac space

$\displaystyle D_{\rm plaq}={\rm
tr_{D}}[(1-\gamma^{\mu})(1-\gamma^{\nu})(1+\gamma^{\mu})(1+\gamma^{\mu})]=-8,$
(57)

where ${\rm tr_{D}}$ means the trace over the Dirac indices. The factor $2$ in
Eq. (56) comes from two directions for each trajectory, which have to be
distinguished in the HPE. The factor $1/n$ in Eq. (27) is canceled by the
number of four starting points of a trajectory; this cancellation occurs for
all trajectories.

Next, let us consider the Wilson loops of length $6$ without windings along
the temporal direction. At this order there are three types of trajectories;
$1\times 2$ rectangle, chair, crown, which are shown in Fig. 2. We define the
operators, $\hat{W}_{\rm rect}$, $\hat{W}_{\rm chair}$, $\hat{W}_{\rm crown}$,
corresponding to these trajectories as

$\displaystyle\hat{W}_{\rm rect}=$ $\displaystyle\frac{1}{N_{\rm c}M_{\rm
rect}N_{\rm site}}\sum_{x}\sum_{\mu\neq\nu}{\rm Re}\ {\rm
tr_{C}}\left[U_{x,\mu}U_{x+\hat{\mu},\mu}U_{x+2\cdot\hat{\mu},\nu}U^{\dagger}_{x+\hat{\mu}+\hat{\nu},\mu}U^{\dagger}_{x+\hat{\nu},\nu}U^{\dagger}_{x,\nu}\right],$
(58) $\displaystyle\hat{W}_{\rm chair}=$ $\displaystyle\frac{1}{N_{\rm
c}M_{\rm chair}N_{\rm
site}}\sum_{x}\sum_{\mu,\nu<\rho,\nu\neq\mu\neq\rho}\sum_{s,t=\pm 1}{\rm Re}\
{\rm
tr_{C}}\left[U_{x,s\nu}U_{x+s\hat{\nu},\mu}U_{x+\hat{\mu},s\nu}^{\dagger}U_{x+\hat{\mu},t\rho}U_{x+t\hat{\rho},\mu}^{\dagger}U_{x,t\rho}^{\dagger}\right],$
(59) $\displaystyle\hat{W}_{\rm crown}=$ $\displaystyle\frac{1}{N_{\rm
c}M_{\rm crown}N_{\rm site}}\sum_{x}\sum_{\mu<\nu<\rho}\sum_{s,t=\pm 1}{\rm
Re}\ {\rm
tr_{C}}\left[U_{x,\mu}U_{x+\hat{\mu},s\nu}U_{x+\hat{\mu}+s\hat{\nu},t\rho}U_{x+s\hat{\nu}+t\hat{\rho},\mu}^{\dagger}U_{x+t\hat{\rho},s\nu}^{\dagger}U_{x,t\rho}^{\dagger}\right],$
(60)

with $U_{x,-\mu}=U^{\dagger}_{x-\hat{\mu},\mu}$ and
$U_{x,-\mu}^{\dagger}=U_{x-\hat{\mu},\mu}$ for $\mu>0$. Eqs. (58)–(60) are
defined so that $\langle\hat{W}_{\rm rect}\rangle=\langle\hat{W}_{\rm
chair}\rangle=\langle\hat{W}_{\rm crown}\rangle=1$ in the weak coupling limit
again; to satisfy this condition, the operators are divided by $M_{\rm
rect}=12$, $M_{\rm chair}=48$, and $M_{\rm crown}=16$, respectively,
corresponding to the number of trajectories per lattice site.

The contribution of these trajectories to the HPE of $\ln\det M(\kappa)$ is
given by

$\displaystyle-2N_{\rm c}\sum_{s}M_{s}D_{s}N_{\rm site}\kappa^{6}\hat{W}_{s},$
(61)

with $s=$rect, chair, crown. $D_{s}$ is the coefficient from the trace in the
Dirac space. For the $1\times 2$ rectangle we have

$\displaystyle D_{\rm rect}=$ $\displaystyle\,{\rm
tr_{D}}\big{[}(1-\gamma_{\mu})(1-\gamma_{\mu})(1-\gamma_{\nu})$
$\displaystyle\times(1+\gamma_{\mu})(1+\gamma_{\mu})(1+\gamma_{\nu})\big{]}=-32,$
(62)

and similar manipulations lead to $D_{\rm chair}=D_{\rm crown}=-16$.

Next we consider the Polyakov-loop type operators having a winding along the
temporal direction. The lowest-order contribution among them is the Polyakov
loop, i.e. the straight lines of length $N_{t}$. To calculate its
contribution, one needs to pay a special attention to the fact that there is
only one independent Polyakov loop for each spatial coordinate on one time
slice, not for each lattice site. Therefore, their contributions to the HPE is
given by

$\displaystyle 2N_{\rm c}D_{\rm pol}N_{s}^{3}\kappa^{N_{t}}\hat{\Omega}_{\rm
R},$ (63)

where the factor $-1$ is to be applied because of the anti-periodic boundary
condition of the quark determinant. The real part of $\hat{\Omega}$ has to be
taken after multiplying the factor 2 because two directions of a trajectory
are independently taken into account. The factor from the Dirac trace for the
Polyakov loop is calculated to be

$\displaystyle D_{\rm pol}={\rm
tr_{D}}\big{[}(1-\gamma_{4})^{N_{t}}\big{]}=2^{N_{t}+1}.$ (64)

Finally, we consider the contribution of the bent Polyakov loops shown in Fig.
3, whose explicit definitions are given by

$\displaystyle\hat{\Omega}_{1}$ $\displaystyle=$ $\displaystyle\frac{1}{N_{\rm
c}M_{\rm bent}N_{\rm site}}\sum_{x}\sum_{i=1}^{3}\sum_{s=\pm 1}{\rm
tr_{C}}\Big{[}\;U_{x,si}U_{x+s\hat{i},4}U_{x+\hat{4},si}^{\dagger}U_{x+\hat{4},4}U_{x+2\cdot\hat{4},4}\cdots
U_{x+(N_{t}-1)\cdot\hat{4},4}\Big{]},$ (65) $\displaystyle\hat{\Omega}_{2}$
$\displaystyle=$ $\displaystyle\frac{1}{N_{\rm c}M_{\rm bent}N_{\rm
site}}\sum_{x}\sum_{i=1}^{3}\sum_{s=\pm 1}{\rm
tr_{C}}\Big{[}\;U_{x,si}U_{x+s\hat{i},4}U_{x+s\hat{i}+\hat{4},4}U_{x+2\cdot\hat{4},si}^{\dagger}U_{x+2\cdot\hat{4},4}\cdots
U_{x+(N_{t}-1)\cdot\hat{4},4}\Big{]},$ (66)

and so on. The factor $M_{\rm bent}=6$ is needed to make
$\langle\hat{\Omega}_{n}\rangle=1$ in the weak coupling limit. From the
definition of $\hat{\Omega}_{n}$ we have
$\hat{\Omega}_{n}=\hat{\Omega}_{N_{t}-n}$. Also, when $N_{t}$ is even
$\hat{\Omega}_{N_{t}/2}$ counts each trajectory twice, and thus its
contribution to the HPE should be divided by $2$. Bearing these facts in mind,
the contribution from $\hat{\Omega}_{n}$ to the HPE of $\ln\det M(\kappa)$ is
given by

$\displaystyle 2N_{\rm c}M_{\rm bent}D_{\rm bent}N_{\rm site}\kappa^{N_{t}+2}$
$\displaystyle\qquad\times\Big{(}\sum_{n=1}^{N_{t}/2-1}{\rm
Re}\,\hat{\Omega}_{n}+\frac{1}{2}{\rm Re}\,\hat{\Omega}_{N_{t}/2}\Big{)}$ (67)

where the last term ${\rm Re}\,\hat{\Omega}_{N_{t}/2}$ should be omitted for
odd $N_{t}$. The contribution of the Dirac trace is calculated to be $D_{\rm
bent}=2^{N_{t}+1}$.

Accumulating these results gives Eqs. (29) and (30).

##  Appendix C Comparison of LO and NLO results

Figure 16:  Binder cumulant $B_{4}$ as a function of $\lambda$ calculated at
the LO. The results at the NLO are also plotted by the thin dashed lines. The
square with errors shows the fit result with three largest volumes. The NLO
result on the same fit result is also shown by the thin circle.

In this appendix, to see the convergence of the HPE at the NLO, we repeat the
analyses in Sec. IV.2 at the LO and compare its results with those at the NLO.
In Fig. 16, we show the Binder cumulant $B_{4}^{\Omega}$ obtained at the LO
along the transition line. In this figure, the NLO results of Fig. 10 are also
shown by the thin dashed lines. We find that, though the difference between
the LO and the NLO results grows as $\lambda$ becomes larger, the deviation is
within a few percent level around the crossing point $\lambda_{\rm c}$.

In Fig. 16, we also show the result of the four parameter fit with Eq. (41)
using data at the LO on $LT=12$, 10 and 9 lattices by the black square. The
same result at the NLO is shown by gray circle for comparison. We find that
the LO result of the fit is consistent with the NLO result within statistical
errors: The values of $b_{4}=1.630(24)$ and $\nu=0.620(47)$ at the LO are
hardly changed from the corresponding NLO values given in Table 3. Though the
central value of $\lambda_{\rm c}=0.00516(15)$ at the LO is about $2.6\%$
larger than the NLO value, it is consistent with the NLO result within errors,
suggesting that the truncation error of the HPE at the NLO is well suppressed
in these quantities.

The success of the $B_{4}^{\Omega}$ fit together with the consistency of the
fit results with the $Z(2)$ values suggests that the $Z(2)$ scaling is
realized also with the LO action when the system volume is sufficiently large.
This is reasonable since the scaling properties near the CP are insensitive to
detailed structures of the theory.

## References

  * Aoki  _et al._ [2006] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, Nature 443, 675 (2006), arXiv:hep-lat/0611014 .
  * Ejiri  _et al._ [2009] S. Ejiri, F. Karsch, E. Laermann, C. Miao, S. Mukherjee, P. Petreczky, C. Schmidt, W. Soeldner, and W. Unger, Phys. Rev. D 80, 094505 (2009), arXiv:0909.5122 [hep-lat] .
  * Asakawa and Yazaki [1989] M. Asakawa and K. Yazaki, Nucl. Phys. A 504, 668 (1989).
  * Stephanov  _et al._ [1999] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. D 60, 114028 (1999), arXiv:hep-ph/9903292 .
  * Asakawa and Kitazawa [2016] M. Asakawa and M. Kitazawa, Prog. Part. Nucl. Phys. 90, 299 (2016), arXiv:1512.05038 [nucl-th] .
  * Bzdak  _et al._ [2020] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, Phys. Rept. 853, 1 (2020), arXiv:1906.00936 [nucl-th] .
  * Bluhm  _et al._ [2020] M. Bluhm _et al._ , Nucl. Phys. A 1003, 122016 (2020), arXiv:2001.08831 [nucl-th] .
  * Pisarski and Wilczek [1984] R. D. Pisarski and F. Wilczek, Phys. Rev. D 29, 338 (1984).
  * Brown  _et al._ [1990] F. R. Brown, F. P. Butler, H. Chen, N. H. Christ, Z.-h. Dong, W. Schaffer, L. I. Unger, and A. Vaccarino, Phys. Rev. Lett. 65, 2491 (1990).
  * Gavin  _et al._ [1994] S. Gavin, A. Gocksch, and R. D. Pisarski, Phys. Rev. D 49, R3079 (1994), arXiv:hep-ph/9311350 .
  * Iwasaki  _et al._ [1996a] Y. Iwasaki, K. Kanaya, S. Kaya, S. Sakai, and T. Yoshie, Z. Phys. C 71, 343 (1996a), arXiv:hep-lat/9505017 .
  * Iwasaki  _et al._ [1996b] Y. Iwasaki, K. Kanaya, S. Kaya, S. Sakai, and T. Yoshie, Phys. Rev. D 54, 7010 (1996b), arXiv:hep-lat/9605030 .
  * Aoki  _et al._ [1999] S. Aoki _et al._ (JLQCD), Nucl. Phys. B Proc. Suppl. 73, 459 (1999), arXiv:hep-lat/9809102 .
  * Karsch  _et al._ [2001] F. Karsch, E. Laermann, and C. Schmidt, Phys. Lett. B 520, 41 (2001), arXiv:hep-lat/0107020 .
  * de Forcrand and Philipsen [2003] P. de Forcrand and O. Philipsen, Nucl. Phys. B 673, 170 (2003), arXiv:hep-lat/0307020 .
  * Bernard  _et al._ [2005] C. Bernard, T. Burch, E. B. Gregory, D. Toussaint, C. E. DeTar, J. Osborn, S. Gottlieb, U. M. Heller, and R. Sugar (MILC), Phys. Rev. D 71, 034504 (2005), arXiv:hep-lat/0405029 .
  * de Forcrand and Philipsen [2007] P. de Forcrand and O. Philipsen, JHEP 01, 077 (2007), arXiv:hep-lat/0607017 .
  * Cheng  _et al._ [2007] M. Cheng _et al._ , Phys. Rev. D 75, 034506 (2007), arXiv:hep-lat/0612001 .
  * Jin  _et al._ [2013] X.-Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda, and A. Ukawa, Phys. Rev. D 88, 094508 (2013), arXiv:1307.7205 [hep-lat] .
  * Ejiri  _et al._ [2016] S. Ejiri, R. Iwami, and N. Yamada, Phys. Rev. D 93, 054506 (2016), arXiv:1511.06126 [hep-lat] .
  * Bazavov  _et al._ [2017] A. Bazavov, H. T. Ding, P. Hegde, F. Karsch, E. Laermann, S. Mukherjee, P. Petreczky, and C. Schmidt, Phys. Rev. D 95, 074505 (2017), arXiv:1701.03548 [hep-lat] .
  * Jin  _et al._ [2017] X.-Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda, and A. Ukawa, Phys. Rev. D 96, 034523 (2017), arXiv:1706.01178 [hep-lat] .
  * Kuramashi  _et al._ [2020] Y. Kuramashi, Y. Nakamura, H. Ohno, and S. Takeda, Phys. Rev. D 101, 054509 (2020), arXiv:2001.04398 [hep-lat] .
  * Banks and Ukawa [1983] T. Banks and A. Ukawa, Nucl. Phys. B 225, 145 (1983).
  * DeGrand and DeTar [1983] T. A. DeGrand and C. E. DeTar, Nucl. Phys. B 225, 590 (1983).
  * Alexandrou  _et al._ [1999] C. Alexandrou, A. Borici, A. Feo, P. de Forcrand, A. Galli, F. Jegerlehner, and T. Takaishi, Phys. Rev. D 60, 034504 (1999), arXiv:hep-lat/9811028 .
  * Saito  _et al._ [2011] H. Saito, S. Ejiri, S. Aoki, T. Hatsuda, K. Kanaya, Y. Maezawa, H. Ohno, and T. Umeda (WHOT-QCD), Phys. Rev. D 84, 054502 (2011), [Erratum: Phys.Rev.D 85, 079902 (2012)], arXiv:1106.0974 [hep-lat] .
  * Fromm  _et al._ [2012] M. Fromm, J. Langelage, S. Lottini, and O. Philipsen, JHEP 01, 042 (2012), arXiv:1111.4953 [hep-lat] .
  * Saito  _et al._ [2014] H. Saito, S. Ejiri, S. Aoki, K. Kanaya, Y. Nakagawa, H. Ohno, K. Okuno, and T. Umeda, Phys. Rev. D 89, 034507 (2014), arXiv:1309.2445 [hep-lat] .
  * Ejiri  _et al._ [2020] S. Ejiri, S. Itagaki, R. Iwami, K. Kanaya, M. Kitazawa, A. Kiyohara, M. Shirogane, and T. Umeda (WHOT-QCD), Phys. Rev. D 101, 054505 (2020), arXiv:1912.10500 [hep-lat] .
  * Cuteri  _et al._ [2021] F. Cuteri, O. Philipsen, A. Schön, and A. Sciarra, Phys. Rev. D 103, 014513 (2021), arXiv:2009.14033 [hep-lat] .
  * Pelissetto and Vicari [2002] A. Pelissetto and E. Vicari, Phys. Rept. 368, 549 (2002), arXiv:cond-mat/0012164 .
  * Ferrenberg and Swendsen [1989] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 63, 1195 (1989).
  * Binder [1981] K. Binder, Z. Phys. B 43, 119 (1981).
  * Rothe [1992] H. J. Rothe,  _Lattice gauge theories: An Introduction_ , Vol. 43 (1992).
  * Iwami  _et al._ [2015] R. Iwami, S. Ejiri, K. Kanaya, Y. Nakagawa, D. Yamamoto, and T. Umeda, Physical Review D 92 (2015), 10.1103/physrevd.92.094507.
  * Hasenfratz  _et al._ [1983] P. Hasenfratz, F. Karsch, and I. O. Stamatescu, Phys. Lett. B 133, 221 (1983).
  * Shirogane  _et al._ [2021] M. Shirogane, S. Ejiri, R. Iwami, K. Kanaya, M. Kitazawa, H. Suzuki, Y. Taniguchi, and T. Umeda (WHOT-QCD), PTEP 2021, 013B08 (2021), arXiv:2011.10292 [hep-lat] .