Thesis for doctoral degree

Canonical approach— Investigation of finite density QCD phase transition —Shotaro Oka
Rikkyo University

Abstract
QCD is the fundamental theory which describes the dynamics of quarks andgluons. If we understand the dynamics at finite density and temperature, i.e. QCD
phase diagram and equation of states, we can progress many studies such as thestudies of unstable nucleus, nuclear fusion, early universe and neutron stars. The
study of QCD phase diagram is very interesting, but we have not understood it well
for a long time. This is because we face a problem in this thesis at finite density. The
problem is called the sign problem. It causes a decrease of the calculation accuracy.
That is why we can not calculate physical quantities with high accuracy at finite
chemical potential.
In this thesis, we try to beat the sign problem using the canonical approachof finite density lattice QCD. Although it is known that the canonical approach
has several numerical problems, we can reduce them and calculate thermodynamic
observables well at finite density. Concretely, in order to reduce the computation
cost of the fermion determinant, we use the winding number expansion, in order toenhance the calculation accuracy of physical quantities, we adopt the multi precisioncalculation to our program based on the canonical approach. In this thesis, wewill see how to improve the canonical approach and a result of thermodynamic
observables which is related with the QCD phase transition at finite density.
Our result shows that we do not observe the peak which represents the confinement–
deconfinement phase transition in baryon number susceptibility. Therefore, we donot see the QCD phase transition yet. However, in this thesis, we find that canonicalapproach can explore the QCD phase diagram beyond µB /T = 3 (µB is the baryon
chemical potential). That is, we explored the QCD phase structure beyond the validity range of Taylor expansion and reweighting method. This opens a bright windowof study of QCD phase diagram at finite density. With our improvement, canonicalapproach has the possibility for investigation of thermodynamic observables at any
chemical potential.

Contents
1 Introduction
1.1 Property of quarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 QCD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Finite temperature and density QCD . . . . . . . . . . . . . . . . . . . . .

1
1
2
5

2 Basic definitions of finite density lattice QCD
8
2.1 Lattice variables and simple lattice actions . . . . . . . . . . . . . . . . . . 8
2.2 Partition function and observables . . . . . . . . . . . . . . . . . . . . . . 10
3 Basic methods of finite density lattice QCD
3.1 Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Improved Monte Carlo methods . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Importance sampling method . . . . . . . . . . . . . . . . . . . . .
3.2.2 Metropolis method . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 Hybrid Monte Carlo method . . . . . . . . . . . . . . . . . . . . .

13
13
13
13
14
15

4 Sign problem
4.1 Sign problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Reweighting technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Taylor expansion method . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Canonical approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16
16
16
17
18

5 Review for study of QCD phase structure with canonical approach20
5.1 Experiment data and canonical partition function . . . . . . . . . . . . . . 20
5.2 Lattice results for QCD phase structure with canonical approach . . . . . 23
6 Calculation methods of fermion determinant
6.1 Reduction formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Hopping parameter expansion . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Winding number expansion . . . . . . . . . . . . . . . . . . . . . . . . . .

27
27
30
30

7 Problems of canonical approach
7.1 Property of discrete Fourier transform . . . . . . . . . . . . . . . . . . . .
7.2 ∗ An instability of DFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 ∗ Reduce these problems . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31
31
33
36

8 Study of finite density phase transition with canonical approach38
8.1 Thermodynamical observables and partition function . . . . . . . . . . . . 38
8.2 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
8.3 Error estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
8.4 ∗ Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
8.4.1 Canonical partition function and analysis of truncation point . . . 41
8.4.2 Pressure, baryon number density and baryon number susceptibility 42

9 Conclusion and discussion

47

A Grassmann integral of fermion action50

B ∗ Algorithm of winding number expansion

51

C Noise method

51

D ∗ Direct method via winding number expansion without canonical ensembles
53
E∗ Algorithm for estimation of cancellation of significant digits59* means the author’s works.

Parts of this thesis have been published in the following journal articles.
• A. Nakamura, S. Oka and Y. Taniguchi, “QCD phase transition at real chemical
potential with canonical approach,” JHEP 02 (2016) 054.
• R. Fukuda, A. Nakamura and S. Oka, “Canonical approach to finite density QCD
with multiple precision computation,” Phys. Rev. D 93, 094508 (2016).
• S. Oka, “Exploring finite density QCD phase transition with canonical approach—Power of multiple precision computation—,” PoS (LATTICE2015) 067 (2016).
• R. Fukuda, A. Nakamura, S. Oka, S. Sakai, A. Suzuki and Y. Taniguchi, “Beating
the sign problem in finite density lattice QCD,” PoS (LATTICE2015) 208 (2016).
• R. Fukuda, A. Nakamura and S. Oka, “Validity range of canonical approach tofinite density QCD,” PoS (LATTICE2015) 167 (2016).

1

Introduction

1.1

Property of quarks

Our world is constructed by the matter such as quarks and leptons. We can predict their
behavior using several gauge theories, that is Weinberg–Salam theory for electroweak interaction and quantum chromodynamics (QCD, Yang–Mills theory for SU(3)) for strong
interaction.
It is well known that, for electroweak interaction, experimental results have been
consistent with predictions based on the perturbative gauge theory with high accuracy.
On the other hand, for strong interaction, we can not quantitively analyze most of the
QCD phase diagram on a temperature–density plane because perturbative methods do
not work in the region.

Figure 1: Schematic QCD phase diagram
However, we have expected that QCD has a rich phase structure at finite temperatureand finite density. This is because qualitative analysis has given the clues which strongly
support the expectation. For example, they suggest that quarks are confined and form
hadronic bound states at a low temperature and a low density. Since our world is
located on this region, our bodies, this paper, the earth and the other planets are madeof hadronic matters.
At high temperature, T > 102 MeV (= 1012 K)[1], quarks are disentangled fromhadrons. This means that quarks can travel almost freely even though they interact
with gluons. This phase is called quark–gluon plasma (QGP). We consider that the1origin of our world was in this phase because our universe began with high temperature.
Therefore, the study of QGP has been one of the most important subjects to reveal theearly universe for ages.
At high baryon density, µB > 103 MeV (= 1018 kg/m3 )[1, 2], quarks can not distinguish the boundary of baryons (see Fig. 2). Therefore, baryons dissolve into a degenerate
Fermi system of quarks; it is called the quark matter. Such an exotic system may exist
in the central core of the neutron stars.

Figure 2: Diﬀerence of hadronic matters and quark matters. R, G and B are the colorsof quarks.

We have mainly obtained a lot of interesting information about the quark matterby experiments, model calculations and lattice QCD simulations. Among these, lattice
QCD simulations play very important roles. This is because it is the first–principles
calculation based on the discretized theory of QCD. In other words, it can reliably
discuss nature of the world managed by strong interaction.

1.2

QCD

QCD describes the dynamics of quarks and gluons. A quark has a color charge which
is one of quantum numbers. It is classified into three states of red, green and blue. Asan electric charge is origin of the electromagnetic interaction, color charge causes thestrong interaction for quarks. Then, the interaction is propagated by gluons.
Let ψ a be a fermion field and Aµ = Aαµ T α be a gauge field. Here, a is a color index;
it describes a type of quarks (a = 1, 2, 3 correspond to red, green and blue, respectively.)T α (α = 1, · · · , 8) are the generators of SU(3). Gluon field Aµ is a matrix which interacts
with fermions and gluons. Therefore, we write it asαα,abAµ = Aab.
µ (x, y, z, t) = Aµ (x, y, z, t)T

2(1)

Then, QCD Lagrangian is written as follows,
1
L = − Gαµν Gαµν + ψ̄ a (iγ µ Dµ − m)ψ a ,4Gαµν = ∂µ Aαν − ∂ν Aαµ + gf αβγ Aβµ Aγν ,(3)− igAαµ T α,ab ψ b ,(4)aDµ ψ = ∂µ ψa(2)where Gαµν is the strength of the gauge field, Dµ is the covariant derivative and g is thecoupling constant. f αβγ are the structure constant of SU(3). f αβγ and T α are written
as follows using the Gell–Mann matrices λα ,
λα
Tα =
,
2
[
]
T α , T β = if αβγ T γ .

(5)
(6)

Here, since gluons have color indices (that is, gluons are not commutative), it caninteract with itself. We can see this from the first term of Eq. (2),()
Gαµν Gαµν = ∂µ Aαν − ∂ν Aαµ + gf αβγ Aβµ Aγν
)((7)× ∂ µ Aα ν − ∂ ν Aα µ + gf αδϵ Aδ µ Aϵ ν= ∂µ Aαν ∂ µ Aα ν − ∂µ Aαν ∂ ν Aα µ + g∂µ Aαν f αδϵ Aδ µ Aϵ ν− ∂ν Aαµ ∂ µ Aα ν + ∂ν Aαµ ∂ ν Aα µ − g∂ν Aαµ f αδϵ Aδ µ Aϵ ν
+ gf αβγ Aβµ Aγν ∂ µ Aα ν − gf αβγ Aβµ Aγν ∂ ν Aα µ
+ g 2 f αβγ Aβµ Aγν f αδϵ Aδ µ Aϵ ν .

(8)The underlined terms of Eq. (8) describe the interactions between gluons, and the otherterms describe the motions of gluons. In the QED, photons do not interact with itself
because the structure constant of U(1) gauge theory satisfies f αβγ = 0.
This characteristic property aﬀects the dynamics of quarks significantly. For example,
quarks (color charges) emit or absorb the gluons as well as electric charges do for thephotons. Then, although photons are emitted radially, gluons are emitted to only one
direction because of gluon–gluon interaction. That is, emitted gluons form a flux tube,
and they are necessarily absorbed by the other color charge. We may consider that sucha flux tube has a binding energy in order to converge the gluons. Because the energy may
become larger with increasing the distance, the flux is stretched between two charges inthe shortest way, and a quark can not appear alone since it needs infinite energy to befree.
In fact, if we try to take away a quark from a baryon, a pair creation of quarks occurs
in the flux tube, we finally get one meson and one new baryon. (See also Fig. 3; it is ameson case.) This phenomenon is called the confinement of quarks.

3

Figure 3: Schematic drawing of gluon flux tube
QCD also has other interesting properties, for instance, the asymptotic freedom is asignificant feature of QCD.
Coupling constant g of QCD is not a constant. In practice, it depends on energy
scale µ which is used for a renormalization, e.g. the momentum of external lines of a
Feynmann diagram. The behavior of g(µ) is given by a calculation of the beta function.
For QCD, It is known that beta function B associates coupling constant g(µ) with energy
scale µ as follows,dg(µ)= −B0 g 3 − B1 g 5 + · · · ,
dµ
(
)
1
11
2
B0 =
Nc − Nf ,(4π)2 33()
1
34 2 10
Nc 2 − 1
B1 =Nc − Nc Nf −Nf ,(4π)2 33Nc

B(g(µ)) = µ

(9)
(10)
(11)

using the two–loop perturbation theory[3]. Here, Nc = 3 is the number of colors and Nf
is the number of flavors. From these, we find the following relation,log(µ) =1
B1
B1+
log(B0 + B1 g 2 ) + (Const.) .
log(g 2 ) −2B0 g 2 2B0 2
2B0

(12)

This equation indicates that µ becomes large when g is small. In other words, decreasing
distance scale a = 1/µ makes coupling constant g small. It means that quarks behave
as asymptotic free particles in a high energy collision although they are confined at low
energy scale. This behavior indicates the presence of QGP at high temperature, andasks us what happens between hadron gas phase and QGP phase; these two phases arerespectively ruled by diﬀerent particles and physics.
The asymptotic freedom was found by a collider experiment of SLAC[4, 5] in 1968,
which was explained theoretically by D.J. Gross and F. Wilczek[6], and H.D. Politzer[7]in 1973. In addition, it is confirmed that the results of string tension*1 of a lattice QCD*1
We can consider the gluon flux tube which bonds two quarks as a string. Then, introducing thestring tension, we can explain the confinement of quarks as a string model. The asymptotic freedom is
explained as the string tension becomes small at short distance scale.

4Figure 4: Running coupling constant dependence of string tension. The left panel andright panel show the results for SU(2) and SU(3), respectively. Vertical axis denotes thequantities χ = a2 σ which is related with the beta function, and horizontal axis denotes
the running coupling constant β = 4/g 2 (left panel) or 6/g 2 (right panel). This figure istaken from Ref. [9].
simulation reproduce the theoretical expectation with high accuracy as Fig. 4. That is,
the non–perturbative result supports the perturbative exception and shows the behaviorof running coupling constant at middle region where the perturbation theory does notwork; it is a valuable result. This simulation was performed by M. Creutz[8, 9] in 1980.

1.3

Finite temperature and density QCD

When we discuss actual QCD systems, we use methods of statistical mechanics andthermodynamics. QCD is a quantum theory, that is, we can apply the path integral
formulation to it. Thus, we can discuss the system of QCD statistically. Thermodynamic
variables, e.g. temperature and chemical potential can be adopted into QCD by thefollowing way.
In order to introduce temperature to QCD, we consider a partition function,()
Z = Tr e−β(Ĥ−µN̂ ) .
(13)
where β = 1/T is an inverse temperature, Ĥ is the Hamiltonian operator and N̂ is thenumber operator. Let |ϕ⟩ be a state of the system on the Hilbert space. Then, we can5rewrite the trace as a summation of physical states,⟩∑⟨Z=
ϕ(x) e−β(Ĥ−µN̂ ) ϕ(x) .

(14)

ϕ

For simplicity, we will neglect fermions and µN̂ for a moment. Then, we notice that theoperator of Eq. (14), exp(−β Ĥ), corresponds to the Euclidean time evolution operatorof quantum mechanics, exp(−Ĥt). We therefore can consider that, essentially, β is theelapsed time of states.
Now, states are represented by gauge field Aµ . Taking the Coulomb gauge, we canwrite the partition function as⟩∑⟨Z=A(β, x) e−β Ĥ A(0, x) ,(15)
A

where 3–dimensional vector A means the spatial components of Aµ . In this equation,
the final state of the bracket is not equivalent to the initial state of it unlike Eq. (14).
In order to remove this inconsistency, we require the periodic boundary condition along
time direction,A(β, x) = A(0, x) .
(16)
Requiring this condition, we can construct the finite temperature QCD. (See also Ref. [10].)In order to introduce density to QCD, we need to consider how to include chemical
potential µ in it. Because we need fermions when we discuss the density of the system,µ should be included in the fermionic part of the Lagrangian, especially in the Dirac
operator.
Then, we can get a hint from Eq. (14). Using the path integral formalism, we canwrite the partition function with chemical potential µ as follows,( ∫ β ∫[])∫
1Z = DADψDψ̄ exp −dτ d3 x − F 2 + ψ̄(γD − m)ψ − µψ † ψ
,
(17)
4
0
where τ = it is the Euclidean time coordinate with the periodic boundary condition.
(Now, we∫ omit the indices.) Here, we use the fact that the number operator is written
as N̂ = d3 xψ † ψ and N̂ is constant of time.
However, it is complicated because ψ † appears. We can simplify this equation via
the change of variablesiAν → Aν − (∂ν eµτ ) e−µτ .
(18)g
The partition function now becomes[( ∫ β ∫])∫1 2
3dτ d x − F + ψ̄(γD − m)ψ
Z = DADψDψ̄ exp −
.
(19)
4
0
That is, the chemical potential appears as the imaginary part of temporal component ofgauge field.
6

A. Roberge and N. Weiss considered above discussion with pure imaginary chemical
potential µ = iµI [11]. Then, because eiµI is periodic, Z(µI ) has the periodicity,()
2πkZ(µI ) = Z µI +,(20)
Nfor any integer k. This result is proved for SU(N) theory. (N = 3 is QCD case.) Thus,
the partition function Z(µ) of QCD has Z3 symmetry which is the center of SU(3)at imaginary chemical potential. This is called Roberge–Weiss (RW) symmetry. In
addition, a phase transition which RW symmetry causes is called RW phase transition.
Using above methods, we can investigate QCD phase diagram at finite temperature and density. Ones may perform it numerically with lattice QCD simulations
and/or several models because it is very hard to investigate QCD theoretically with
non–parturbative methods.
However, we have not been able to study it well, although we use the numericalmethods. For example, the QCD critical end point have been searched for a long time.
It is known that, in T –µ plane, the critical point is located at Tc ≈ 150 MeV from zerodensity lattice QCD simulations*2 [1]. But, critical density (chemical potential) µcB is
not known well.
Z. Fodor and S.D. Katz studied it in 2004 for 2 + 1 flavor lattice QCD[12], then theygave µcB = 320(40) MeV. On the other hand, S. Datta, R.V. Gavai and S. Gupta studied
it in 2012 for two flavor lattice QCD[13]. They then gave µcB = 243(20) MeV. Moreover,
O. Scavenius, A. Mocsy, I.N. Mishustin and D.H. Rischke introduced the phase diagramsof σ–model and NJL model which are eﬀective models of QCD, in 2000[14]. From thispaper, critical point of these models are at µcB = 220 MeV and 330 MeV, respectively.
As we saw, the location of the critical point depends on models. In addition, at a
finite density, there is the sign problem which is discussed in section below. This problem
causes the decease of calculation accuracy. That is why lattice QCD simulations (andseveral models) can not estimate physical quantities well in this region of finite chemical
potential.
Although the sign problem is a severe diﬃculty,ones have been studying how to avoid
this problem. As one of such methods, the canonical approach was developed. Thismethod has potential to avoid the problem, but the other several numerical problems
are left in the method.
In this thesis, we try to completely avoid the eﬀect of the sign problem and investigate
the QCD phase diagram at finite density by solving the problems of canonical approach.
*2

The critical point is located at finite density. This means that we can not use usual order parameterssuch as the Polyakov loop and chiral condensate because the center symmetry of SU(3) and chiralsymmetry is broken by the massive fermions. Thus, when we discuss the location of critical point, we
use the critical temperature at zero density alternatively or pseudo critical temperature Tpc which isapproximately estimated at finite density by the Polyakov loop or chiral condensate. In this thesis, wewrite the both of critical and pseudo critical temperature as Tc .

7

2

Basic definitions of finite density lattice QCD

2.1

Lattice variables and simple lattice actions

Let us begin to consider the lattice QCD. In this theory, space–time coordinates are
discretized as lattice. We put fermions and gluons on the lattice. More precisely, fermionsare on the sites of the lattice, and gluons are on the links of the lattice because they
have a direction as Lorentz index.

Figure 5: Schematic sketch of lattice variables
Let a be the lattice spacing. Because Aµ is not covariant under the SU(3) gaugetransformation, we alternatively use the following covariant quantity,Uµ (x) = U (x → x + µ̂a) := eiagAµ (x+µ̂/2) ∈ SU(3) ,(21)

for a lattice variable. Here, µ̂ is the unit vector along µ direction. Uµ are called the linkvariables.
Since they are unitary matrices, we find the following relation,Uµ (x)U−µ (x + µ̂a) = 1 = Uµ (x)Uµ† (x) ,(22)thus U−µ (x + µ̂a) = Uµ† (x).
We can write the gauge action by gauge invariant quantities which are constructedby Uµ . For example, the trace of closed loops of Uµ is gauge invariant. The simplest
loop, that is a unit square on the lattice,
trUµ (x)Uν (x + µ̂)Uµ† (x + ν̂)Uν† (x) ,8(23)Figure 6: Schematic drawing of plaquetteis called the plaquette. It is invariant under the following gauge transformation,
Ψ(x) → Ψ′ (x) = e−iθ (x)T Ψ(x) =: Ω(x)Ψ(x) ,iAµ (x) → A′µ (x) = Ω(x)Aαµ T α Ω† (x) − [∂µ Ω(x)]Ω† (x) ,
g
α

α

(24)
(25)

with Ω(x) ∈ SU(3). That is,
Uµ (x) → Uµ′ (x) = Ω(x)Uµ (x)Ω† (x + µ̂a) .

(26)

Using the plaquette, we can write the gauge action as
SG =β∑trUµ (x)Uν (x + µ̂)Uµ† (x + ν̂)Uν† (x) .
3

(27)

µ̸=ν

Now, we use β as the eﬀective coupling constant. Note that it is not the inverse temperature. This action is called Wilson gauge action.
It is diﬃcult to construct the fermion action. This is because fermions get extra
degrees of freedom when we discretize the Dirac operator naively. This problem isknown as the fermion doubling problem.
To reduce the problem, K.G. Wilson considered the following action,∑
SF = a4
Ψ̄(x)∆(x, x′ )Ψ(x′ ) ,
(28)
x,x′

∆(x, x′ ) = δx,x′ − κ

4 (
∑
(

)
)
(
)
r − γi Ui (x)δx′ ,x+î + r + γi Ui† (x′ )δx′ ,x−î ,

i=1

9

(29)

where κ is the hopping parameter. κ rules the movement of quarks, which depends on
the quark mass mq [15] as
1
κ=
.
(30)
8 + 2mq a
The parameter r ranges from zero to one. Nonzero r gives the mass which diverges at
a = 0 to the extra fermions. We usually set r = 1. This term is called Wilson term.
Although Wilson terms break the chiral symmetry on the lattice because they introduce
extra masses into the theory, the symmetry is recovered when we take the continuum
limit a → 0. Such fermion action is called Wilson fermion action.

2.2

Partition function and observables

In order to calculate physical quantities, let us firstly consider a representation of the
grand partition function. We can construct it using path integral formulation,∫
Z(µ) = DU DΨDΨ̄e−SG e−SF (µ) .
(31)
Then U is link variables, Ψ and Ψ̄ are fermions and anti–fermions on a lattice respectively.
SG and SF are a gauge action and a fermion action respectively. We will see them indetail below.
For the gauge part SG , we use the improved Wilson gauge action,
SG =

]
β ∑[
1×1
1×2
c0 trPµν+ c1 trPµν
,
3

(32)

µ̸=ν

1×1
where β = 6/g 2 is the eﬀective coupling constant (g is a bare coupling constant), Pµν
1×2
is the usual plaquette and Pµν is the rectangle plaquette:
1×1Pµν(x) = Uµ (x)Uν (x + µ̂)Uµ† (x + ν̂)Uν† (x)

(33)

1×2
Pµν
(x) = Uµ (x)Uµ (x + µ̂)Uν (x + 2µ̂)Uµ† (x + µ̂ + ν̂)Uµ† (x + ν̂)Uν† (x) .

(34)

We set the coeﬃcients in Eq. (32) as c1 = −0.331 and c0 = 1 − 8c1 as Ref. [16]. It iscalled Iwasaki gauge action.
This action was suggested by Y. Iwasaki in 1983. The aim of this action is that thenumerical errors from lattice spacing a are reduced when we calculate the action. In
other word, the lattice action approaches the original continuous action rapidly when we
take the continuum limit a → 0.
The fermion part SF is written as∑Ψ̄(x)∆(µ, x, x′ )Ψ(x′ )SF (µ) = a4
x,x′

10

(35)

Figure 7: Plaquette and rectangle plaquette. The left panel and right panel is the usual(unit square) plaquette and rectangle plaquette, respectively.
where ∆(µ) is a ferimon matrix. Now, we adopt the improved Wilson fermions withchemical potential µ for future reference,′∆(µ, x, x ) = δx,x′ − κ

3 (
∑
(

)
)
(
)
1 − γi Ui (x)δx′ ,x+î + 1 + γi Ui† (x′ )δx′ ,x−î

i=1

(
)
(
)
(
)
− κ e+µa 1 − γ4 U4 (x)δx′ ,x+4̂ + e−µa 1 + γ4 U4† (x′ )δx′ ,x−4̂∑− κCSW δx,x′
σµν Fµν (x) ,

(36)

µ<ν

where the second line which includes e±µa denotes the quarks hopping along time direction (that is, 4̂ represents the unit vector along time direction),
1×1
1×1
1×1
1×1
Fµν = Pµν(x) + Pν−µ(x) + P−µ−ν(x) + P−νµ
(x) ,

(37)

and σµν = i[γµ , γν ]/2.
Then, we introduce eµa to the terms of U4 and U4† because of Eq. (18). This meansthat fermions get the coeﬃcient eµa when they travel along t–coordinate. Since theinverse temperature is related with the boundary conditions of ψ and Aµ , if fermionswind around the lattice along the time direction, they get eµaNt = eµ/T .
The term of CSW is called the clover term. This is the correction term of Wilson
fermions, which includes the field strength on the lattice Fµν . The term reduces thediscretization error of the action. By a result of one–loop perturbation theory[17], weset CSW = (1 − 0.8412/β)−3/4 . (This β is the coupling constant.)
Note that we can integrate out e−SF because it forms the Gaussian integral of theGrassmann number Ψ. (See Appx. A.) As a result, Eq. (31) is rewritten as follows.
∫
Z(µ) = DU [det ∆(µ)] e−SG .
(38)

11

Figure 8: Schematic drawing of clover term.
It is obtained for one flavor fermion case. If we want to consider many flavor fermions,
we can use the following equation,∫(39)
Z(µ) = DU [det ∆(µ)]Nf e−SG ,where Nf is the number of fermions.
Using above definitions, we can calculate the expectation values of thermodynamic
observables,∫1⟨O(µ)⟩ =DU O(µ)[det ∆(µ)]Nf e−SG .
(40)
Z(µ)
Here, operater O(µ), fermion matrix ∆(µ) and gauge action SG depend on link variables
U . Therefore, in a practical caluculation, we calculate them for each configuration (aset of {U }) and average their product.

12

3

Basic methods of finite density lattice QCD

3.1

Monte Carlo method

Here, link variables U have many indices, specificallyU = Uµab (x, y, z, t) ,(41)

where a and b are the color indices, µ is the Lorentz index, and (x, y, z, t) are the space–time indices.
”DU ” means the integration with all of the combinations of U . Now, a, b = 1, 2, 3
and µ = 1, 2, 3, 4, and we discretize the coordinates x, y, z, t = 1, · · · , 10 for instance.
Then, DU represents 3.6 × 105 times∫ integration.
Moreover, we have to perform DUµab (x, y, z, t) using the division quadrature. Thatis, if we write the component of gauge field with fixed a, b, µ and (x, y, z, t) as Uµab (x, y, z, t) =:
u,
∫ (
∫
)
d Uµab (x, y, z, t) f [U ] = duf (u)≈N∑

f (uk )∆u ,

(42)

k=1

where uk (k ∈ {1, 2, · · · N }) is a sampled u. Therefore, if N = 10, we should perform
5
103.6×10 operations in this case. Hence, it is diﬃcult to calculate Eq. (40) naively.
We then can use the Monte Carlo method in order to reduce the computational costof the integration. In this method, we do not use the sequence number, and alternatively
use the random number as U . Here, we consider the probability density function p of
the random number asp ∝ [det ∆(µ)]Nf e−SG ,(43)
in this case. We are going to see it in detail below.

3.2

Improved Monte Carlo methods

3.2.1

Importance sampling method

When we generate the gauge configurations using the Monte Carlo method, we naively
use the uniform random numbers. However, sometimes it is not an eﬃcient method
because it generates not only statistically important configurations but also unimportantothers with the same probability. Therefore, we have to improve it in order to pick only
important ones from them.
Importance sampling method is one of the methods which enhances the statistic
accuracy. In this method, we adopt alternative probability density function and use it13

as a filter. For instance, we consider a simple Monte Carlo case,∫
E[O] = O(x)f (x)dx ,(44)

where E[O] is the expectation value of a function O(x), f (x) is the probability density
function, x is the integration variable. We calculate this in practice,
1 ∑O(xm )
E[O] = lim
m→∞ m
m

(45)

j=1

where xm are the random numbers which conforms to the distribution f (x). Then, if
we want to change the probability density f (x) to another convenient function f ′ (x), we
can rewrite Eq. (44) as∫
f (x) ′
E[O] = O(x) ′f (x)dx .
(46)
f (x)
Thus, Eq. (45) can be rewritten,
1 ∑f (x′ )
E[O] = limO(x′m ) ′ m,
m→∞ mf (x′m )
m

(47)

j=1

where x′m are the random numbers which conforms to the distribution f ′ (x), f (x)/f ′ (x)plays the role of the weight function.
Using such a procedure, if f (x) is intricate or it is not useful function, we can use
any useful functions through f ′ (x). Naively, we use the uniform distribution as f ′ (x),but we usually use more suitable function, e.g. Gauss distribution in order to enhance
the accuracy of the integration.
3.2.2

Metropolis method

To perform the Monte Carlo integration, we have to generate the set of integration
variable x which conforms to the distribution f (x) or f ′ (x). In this thesis, we use theMetropolis method for this purpose.
Let Q(x|y) be a conditional probability density function. It generates new sample x from previous sample y. Q(x|y) is also called the proposal density or jumping
distribution. Now, we require that Q has to satisfy the detailed balance condition,
Q(x|y) = Q(y|x) .

(48)

Metropolis algorithm is constructed as follows.
1. Initialization: choose an arbitrary point x0 as an initial sample and an arbitrary
probability density Q(x|y).
2. Iteration: for each t,
14

(a) Generate new sample x′ from Q(x′ |xt ).
′(x )(b) Calculate acceptance ratio α = ff (x. Note that it is just the weight function.
t)(c) If α ≥ 1, x′ is more likely than xt , then it is accepted for next sample xt+1 ;if α < 1, x′ is accepted according to probability α; when x′ is rejected, set
xt+1 = xt .
We can generate the gauge configurations using this algorithm, but it has a demerit— the generated samples are correlated. This means that, when we use these samples,
we have to extract every n–th samples and only use them (n is determined by the
autocorrelation between adjacent samples. we do not discuss it in detail.) Thus, thisalgorithm is not eﬃcient. Usually, an improved method is alternatively used; we will see
it in the following section.
3.2.3

Hybrid Monte Carlo method

Hybrid Monte Carlo (HMC) method is one of the Markov chain Monte Carlo (MCMC)
methods. It is constructed by the combination of hybrid molecular dynamics (HMD)
and the Metropolis method, and mainly used in the two flavor lattice QCD simulations.
1. Initialization:(a) Choose ϕ as the integration variable.
(b) Generate the conjugation momentum π with probability P (π) ∝ e−π /2 (Gauss
distribution).
2

2. Iterations (HMD part): evolve ϕ and π through the following equations,1 2 H[ϕ, π] = − π + S(ϕ)
2∂H[ϕ, π]
ϕ̇ =
=π

∂π



 π̇ = − ∂H[ϕ, π] = − ∂S(ϕ)
∂ϕ
∂ϕ

(49)

where H[ϕ, π] is the Hamiltonian and S(ϕ) is the action for ϕ (it is given).
3. Metropolis test: using the Metropolis method, accept or reject ϕ and π according
to the following probability,{}
P ({ϕinit , πinit } → {ϕnew , πnew }) = min 1, e−∆H ,(50)∆H = H(ϕnew , πnew ) − H(ϕinit , πinit ) .
4. Repeat 2 and 3.

15

(51)

4

Sign problem

4.1

Sign problem

Basically, using above definitions and methods, we can start the calculation of finite
density lattice QCD. However, det ∆ becomes a complex number when µ is nonzero. We
can see this using the definition of ∆ and the property of the determinant. First, weobserve the following relation by a property of a determinant:[det ∆(µ)]⋆ = det ∆† (µ) .

(52)

Then, ∆† appears, and we can calculate it using Eq. (36). Here, we use the following
equations,(1 − γi )† = 1 − γi(i = 1, 2, 3)= γ5 (1 + γi ) γ5eµ (1 − γ4 )† = e

µ⋆

(53)

(1 − γ4 )µ⋆= e γ5 (1 + γ4 ) γ5 .

(54)Note that we now use the Euclidean gamma matrix {γi , γj } = 2δij . Thus, we find[det ∆(µ)]⋆ = det (γ5 ∆(−µ⋆ )γ5 )(55)⋆= det (γ5 ) det (∆(−µ )) det (γ5 )⋆= det ∆(−µ ) .

(56)
(57)

det ∆(µ) is a complex number if Re(µ) ̸= 0.
When we calculate thermodynamic observables in the finite density region (Re(µ) ̸=0), Monte Carlo method does not work well because the probability [det ∆(µ)]Nf e−SG iscomplex. Such a problem is called the sign problem[18, 19]. This problem was found in
1985[20, 21]. Although researchers have studied it about 30 years, it is not solved yet.

4.2

Reweighting technique

The reweighting technique[22] may be the most famous method which is made in order
to avoid the sign problem. In this method, we rewrite fermion matrix as
det ∆(µ) =det ∆(µ)
det ∆(µ0 ) .
det ∆(µ0 )We usually set µ0 = 0. Then, partition function becomes[]∫det ∆(µ) Nf
Z(µ) = DU[det ∆(µ0 )]Nf e−SGdet ∆(µ0 )][∫1det ∆(µ) Nf= Z(µ0 ) ·[det ∆(µ0 )]Nf e−SG
DU
Z(µ0 )det ∆(µ0 )⟨[] ⟩det ∆(µ) Nf,= Z(µ0 ) ·det ∆(µ0 )
µ0

16

(58)

(59)

where ⟨·⟩µ0 means the expectation value on the gauge configurations which are generated
at µ0 . Similarly, the observables is written as[]
det ∆(µ) Nf
DU O(µ)[det ∆(µ0 )]Nf e−SGdet ∆(µ0 )
[
]
∫
1
1
det ∆(µ) Nf=
Z(µ0 ) ·DU O(µ)[det ∆(µ0 )]Nf e−SGZ(µ)
Z(µ0 )det ∆(µ0 )⟨] ⟩ /⟨[] ⟩[det ∆(µ) Nf
det ∆(µ) Nf= O(µ).
(60)
det ∆(µ0 )det ∆(µ0 )1⟨O(µ)⟩µ =
Z(µ)

∫

µ0

µ0

At µ0 = 0 or pure imaginary, we can calculate ⟨·⟩µ0 because the determinant of fermion
matrix is real. In other words, we consider the complex probability of the Monte Carlointegration as a part of the observables.
Note that this method is just one of the importance sampling method.
This method seems to work out, but it does not work well in practice. The reason is
as follows. When we vary µ, the distribution of det ∆(µ) is varied drastically. Therefore
the distributions do not overlap and we can not sample the gauge configuration well.
This problem is called the overlap problem.
Multi–parameter reweighting (MPR) technique is the expansion method of original
reweighting technique which is made in order to reduce the eﬀect of overlap problem. In
this method, we reweight e−SG too. Concretely,
Nf −SG (β)[det ∆(µ)]

e[]det ∆(µ) Nf e−SG (β)=· [det ∆(µ0 )]Nf e−SG (β0 ) ,det ∆(µ0 )e−SG (β0 )(61)where β is the eﬀective coupling constant (see Eq. (32), the definition of SG ). Using thistechnique, we can enhance accuracy in part at high temperature[23], T > Tc , but theoverlap problem still remains in the other part[24].

4.3

Taylor expansion method

We calculate the thermodynamic observables via the partition function Z(µ). For example, the pressure is given as∆p(µ, T )p(µ, T ) − p(0, T )
:=
4
T
T4
1
[log(Z(µ, T )) − log(Z(0, T ))]=V T3Z(µ, T )
Nt 3
log
=
.
Nx Ny Nz
Z(0, T )17(62)where p is the pressure and V = Nx a × Ny a × Nz a is the spatial volume of the system.
Then, we can expand it in powers of µ/T around µ = 0,∆p(µ, T )
=
T4

∞
∑

cn (T )

n=2,4,6,···( µ )nT.

(63)

Here, the even powers of µ/T only survive because the QCD partition function has thesymmetry Z(µ) = Z(−µ). We can calculate cn via Z(µ) as follows,
cn =
=

1
∂n
∆p(µ, T )nn! ∂(µ/T )
T4
µ=0
1
Nt 3
∂n
.
log Z(µ, T )n! Nx Ny Nz ∂(µ/T )n
µ=0

(64)

Because cn does not depend on µ, we can calculate ∆p/T 4 at any µ using Eq. (63) once
we calculate cn .
In order to use such a method, we have to calculate Z(µ) at several values of µ.
Although it is diﬃcult to perform the calculation at finite µ because of the sign problem,
it is easy to perform near µ = 0. Thus, we calculate cn around µ = 0 and extrapolate
the observables to finite µ. Such method is called Taylor expansion method, which is
employed in order to calculate thermodynamic observables at any µ.
However, in general, although this method works out well at µ/T ≤ 1, it may notwork well at µ/T > 1 because the numerical errors become severe. In fact, the studies of
Taylor expansion method[13, 25, 26, 27] show the results of equation of state for µ/T ≤ 1.
(Several studies[28] show it for µ/T ≤ 1.2.) Here, K. Nagata and A. Nakamura notethat Taylor expansion and MPR methods work out at the same region of µ[24]. They
concluded that they are consistent methods and the overlap problem of MPR and thetruncation error of Taylor expansion method are negligible. That is, we may not explore
QCD phase diagram beyond µ/T = 1 using these methods.
We need to develop another method in order to explore the QCD phase diagram.

4.4

Canonical approach

Canonical approach is one method to avoid the sign problem. It was suggested by A.
Hasenfratz and D. Toussaint[29] in 1992.
They focused on the relation between grand partition function Z(µ) and canonicalpartition function ZC (n). In the statistical mechanics, it is well known that we can write18Z(µ) as a polynomial of fugacity,()
Z(µ) = Tr e−(Ĥ−µN̂ )/T
∞ ⟨
∑

=

n=−∞
∞
∑

=:

⟩
n e−Ĥ/T n enµ/T()n
ZC (n) eµ/T
,

(65)

n=−∞

with a thermodynamical limit. Here, |n⟩ is the eigenvectors of the number operator N̂
(N̂ |n⟩ = n|n⟩) and ZC (n) is the canonical partition function with net quark number n.
In the second line, we use the fact that the number operator does not depend on time,that is [H, N ] = 0. Equation (65) is called the fugacity expansion.
Hasenfratz and Toussaint noticed that Z(µ) is just the Laplace transform of ZC (n).
They consider that, if we use pure imaginary chemical potential µ = iµI (µI is real), we
get Z(µ) as the inverse Fourier transform of ZC (n),Z(µ = iµI ) =∞∑ZC (n)ei(µI /T )n .

(66)

n=−∞

Therefore, we also get ZC (n) as the Fourier transform of Z(µ),
∫ ( )
1
µI
ZC (n) =dZ(µ = iµI )e−i(µI /T )n .
2π
T

(67)

Note that, for pure imaginary chemical potential, det ∆ is real and we can calculateZ(iµI ) and ZC (n) without the sign problem. Moreover, if we get ZC (n), we can calculate
Z(µ) with ANY chemical potential using Eq. (65). Thus, adopting such a procedure,
we can calculate Z(µ) without the sign problem. This method is called the canonicalapproach.
In this thesis, we use the canonical approach in order to avoid the sign problem.

19

5

Review for study of QCD phase structure with canonicalapproach

5.1

Experiment data and canonical partition function

The canonical partition function ZC (n) is related with the number of events of thecollider experiments. This is because ZC (n) represents the existence probability of thesystem which n particles are in. That is, we can get the following equation from Eq. (65),1=∑ ZC (n)(eµ/T )n
Z(µ)

n

.

(68)

From this equation, we can consider ZC (n)(eµ/T )n /Z(µ) as the existence probability.
Here, let us consider the ratio ZC (n)/ZC (0). This ratio is the normalized canonicalpartition function. We will use it when we discuss the existence probability because
(
)0
(
)n /
ZC (0) eµ/T
ZC (n) eµ/T
ZC (n) ( µ/T )n
=
e
.
(69)
Z(µ)
Z(µ)
ZC (0)
Thus, we can cancel Z(µ) using the ratio.
()n
We can see normalized existence probability ZC (n) eµ/T /ZC (0) in collider experiments. For example, Fig. 9 shows the net proton dependence of the number of eventsin the Au+Au collision at RHIC.
Let Nevent (np ) is the number of event with the net proton number np = 3n and Ntotal
is the total number of events. Then, normalized number of events Nevent (np )/Nevent (0)is given as follows,( µ/T )np /
Z(µ)
N
·
Z
(n
)
e
p
C
total
Nevent (np )/=()0Nevent (0)N· Z (0) eµ/T
Z(µ)
C

total

=

Z (n ) (Cp

ZC (0)

eµ/T

)np

.

(70)

Note that ZC (n) satisfies
ZC (n) = ZC (−n) .

(71)

This is because the QCD partition function has the symmetry of the exchange of perticles
for anti–particles (C–symmetry).
Using Eqs. (70) and (71), we can calculate ZC (n)/ZC (0) from RHIC data. Figure 10
shows the net proton number dependence of normalized canonical partition function
from RHIC data. From this figure, we see that ZC (n) becomes smaller rapidly with√increasing |n| and it is raised with increasing the collision energy sN N .

20

Number of Events

10 6

Au+Au 200 GeV

10 5

0.4<pT <0.8 (GeV/c)

0-5%30-40%
70-80%|y|<0.5

10 4
10 3
10 2
10
1
-20

-10

0

10

20

Net Proton ∆(N p)
Figure 9: Net proton dependence of number of events in Au+Au collision at RHIC. This
figure is taken from Ref. [30].

21

1
-4

10

-8

10

(
-12

-10

0

(

10

-30

(

(
(
(

-20

19.6GeV)19.6GeV)
39 GeV) 39 GeV)
200GeV) 200GeV)

n

10

20

30

Figure 10: Proton number dependence of normalization canonical partition function√
from RHIC data with sN N = 19.6, 39 and 200 GeV. Vertical axis represents normalized canonical partition function ZC (np )/ZC (0) and horizon axis represents net proton
number nq . This figure is taken from Ref. [31].

22

5.2

Lattice results for QCD phase structure with canonical approach

In this section, we focus on the lattice QCD results with canonical approach at finite
density. We will not introduce the results for SU(2) or results at zero density.
Because the canonical approach requires large computation resources, we could not
perform full QCD simulation with canonical approach well in 1990’s. In the first study
with canonical approach[29], Hasenfratz and Toussaint performed the calculation of number density on 24 and 44 lattice. See Fig. 11. Their results tell us that, at low temperature, we can not perform this calculation well, although at high temperature, we may
perform it well.

Figure 11: Number density as a function of chemical potential by Hasenfratz and Toussaint in Ref. [29]. The left and right panel is calculated at high and low temperature,
respectively.

In 2004, Y. Sasai, A. Nakamura and T. Takaishi investigated the crossover phasetransition of QCD at finite density[32]. In this study, they used 83 × 4 lattice. They
calculated the chiral condensate as Fig. 12. From this figure, we see the crossover phasetransition as the jump of the chiral condensate. Thus, we can conclude that the canonicalapproach can be used in the study of phase transition. Note that, if Tc ≈ 200 MeV, theystudied at µ ≈ 100 MeV.
Ph. de Forcrand and S. Kratochvila studied the QCD phase diagram via baryon
number density with canonical approach in 2005[33, 34]. They then used 63 × 4 lattice.
Their results are shown in Figs. 13 and 14. In this study, they calculated the baryon
number beyond µ/T = 1. It seems that, using the canonical approach, we can explore
the QCD phase diagram at high density. However, they note that their results areaﬀected by the finite size eﬀect strongly. Therefore, we have to use larger lattice in order
to investigate the QCD phase transition.

23

Figure 12: Chiral condensate as a function of the eﬀective coupling constant at several
finite density by Sasai, Nakamura and Takaishi in Ref. [32].

Figure 13: Baryon number as a function of µ/T by de Forcrand and Kratochvila inRef. [33]. Blue points correspond to their results.

24

Figure 14: Phase diagram in the T –µ plane by de Forcrand and Kratochvila in Ref. [33].
Red points correspond to their results.

25

Figure 15: Comparison between the number density from canonical approach and Taylor
expansion method. This figure was made by Danzer and Gattringer in Ref. [36].
The lattice QCD group of Kentucky (χQCD collaboration) performed several important studies for canonical approach. For example, they developed the winding number
expansion for the fermion matrix in Ref. [35] in 2008. This is introduced below in detail.
Using this method, we can perform the calculation of the fermion determinant fasterthan the other methods.
J. Danzer and C. Gattringer performed the calculation of number density and number
susceptibility at finite chemical potential with 83 × 4 lattice in 2011[36]. Then, they
compare the results of canonical approach with the results of Taylor expansion method
as Fig. 15 and Tab. 1. From Fig. 15, we see that the results of canonical approach andTaylor expansion method are consistent. However, from Tab. 1, we see the diﬀerence
between two methods.
Danzer and Gattringer noted that, at finite volume, fugacity series is a finite Laurent
series in the fugacity parameter eµ/T , while the Taylor series is an infinite series evenat this volume. Since the fugacity series and Taylor series have diﬀerent properties,
they suggested that we have to choose a better method (fugacity expansion or Taylor
expansion) for our purposes.

26

Table 1: Comparison between the number susceptibility from canonical approach (fugacity expansion) and Taylor expansion method. This table was made by Danzer andGattringer in Ref. [36].

6Calculation methods of fermion determinant

Before discussing the problems of canonical approach, we discuss the computing cost of
the fermion determinant.
Since the fermion matrix ∆(µ) has many degrees of freedom, its rank becomes verylarge. Concretely, rank ∆ = Nx Ny Nz Nt × 4 × 3 (Nx is a number of points of x coodinate,
and so on.) and if we choose Nx = Ny = Nz = Nt = 10, rank ∆ ≈ 105 . Moreover, wehave to calculate it for EACH configuration.

6.1

Reduction formula

We firstly introduce the reduction formula of Wilson fermion matrix. It was suggested by
Nagata and Nakamura[37] in 2010. (See also Ref.[38].) Note that the reduction formulafor staggered fermion have been suggested since 1986 [29, 39]. Nagata and Nakamura
noticed that Wilson fermions form a sparse band matrix. It can be resolved into the
diagonal part and its neighbor parts.
For future reference, we divide the determinant of Wilson fermions into three terms,∆ = B − 2z −1 κr− V − 2zκr+ V † .

27

(72)

Here, r± = (1 ± γ4 )/2 is the projection operators, z = e−µa andB(x, x′ ) := δx,x′ − κ

3 (
∑
(
i=1

− κCSW δx,x′))()
1 − γi Ui (x)δx′ ,x+î + 1 + γi Ui† (x′ )δx′ ,x−î∑σµν Fµν (x) ,

(73)

µ<ν
′

V (x, x ) := U4 (x)δx′ ,x+4̂ ,(74)(x, x ) := U4† (x)δx′ ,x−4̂ .

(75)

V

†

′

These matrices are also given as follows in a time–plane block matrix form,
′

t=1
t=2
t=3
..
.

B=

t = Ntt′ = N
t
0
0
0
0 

0
0 
 ,
..
.. 
.
. 

···
t = 1B1 00 ··· 0 B2 0 · · · 0
0 B3 · · ·

 ..
..
..
..
 .
.
.
.
0
0
0 ···

U4 (t = 1)
0
···

0
U4 (t = 2) · · ·


..
..
..
.
.
.
V =


0
0
·
·
·
·
·
·


0
0
···
···
U4 (t = Nt )
0
···
···
0
0
..
.

(76)

0 BNt
0
0
..
.

0
0
..
.







 .

U4 (t = Nt − 2)00U4 (t = Nt − 1) 
0
0

(77)

In order to cancel V † , we multiply ∆ andAs a result,

P := r− + z −1 r+ V .

(78)∆P = (Br− − 2κr+ ) + (Br+ − 2κr− )z −1 V ,(79)Here, (Br− − 2κr+ ) is the diagonal part and (Br+ − 2κr− )z −1 V is the other part. As
the time–plane block matrix, it is given as
α1β1 z −1
α2β2 z −1




..
..
(80)
∆P = 
 ,
.
.


−1


αNt −1 βNt −1 z
−1
βNt z
αNt
where,
αi = αab,µν (⃗x, ⃗y , ti )µν abσν= B ab,µσ (⃗x, ⃗y , ti )r−− 2κr+δ δ(⃗x − ⃗y ) ,βi = β

ab,µν

=B

(81)

(⃗x, ⃗y , ti )

ac,µσ

µν
σν cb
(⃗x, ⃗y , ti )r+
U4 (⃗y , ti ) − 2κr−
δ(⃗x − ⃗y )U4cb (⃗y , ti ) .

28

(82)

Here, a, b is color indices
Using these equations, we can write the determinant of ∆P ,
(N
)
∏t
(
)
det(∆P ) =det(αi ) det 1 + z −Nt Q ,

(83)

i=1
−1
where Q := (α1−1 β1 ) · · · (αNβ ). Since det(P ) = z −N/2 and N = Nx Ny Nz Nt × 3 × 4,
t Nt
we get
(N
)
∏t
(
)
N/2
det ∆ = zdet(αi ) det 1 + z −Nt Q
i=1= z −N/2

(N
∏t

)

(
)
det(αi ) det z Nt + Q .

(84)

i=1

Then, we use a property of determinant det(AB) = det(A) det(B), and det(cA) =cN det(A) with a scalar c and the rank N .
We are going to see that this method is related with the canonical partition function.
Let λn be the eigenvalues of Q. Then, the determinant of the reduced matrix iswritten as
N
red
∏
Nt
(85)
(λn + z Nt ) ,det(z + Q) =n=1

where Nred is the rank of the reduced matrices. We can expand it as
z

−N/2

N
red
∏

∑

Nred /2(λn + z

Nt

)=

n=1

( )n
cn z Nt .

(86)

n=−Nred /2

Here, using z Nt = eµaNt = eµ/T , we finally obtain the following relation,∑

Nred /2det ∆(µ) =()n.
Cn eµ/T(87)

n=−Nred /2

We can calculate Cn numerically. Nagata and Nakamura noticed later that Cn is just the
canonical partition function ZC (n). Thus, we can calculate ZC (n) with the eigenvaluesof reduction matrix.
As seen above, we can dissolve the determinant of ∆ into one of αi and βi . Then,
we reduce the rank of matrix, from rank∆ = Nx Ny Nz Nt × 4 × 3 to rankαi = rankβi =Nx Ny Nz × 4 × 3 = (rank∆)/Nt . The time complexity of a determinant calculation3 ) in a typical algorithm. In the reduction formula,depending on the rank Nrank is O(Nrank
we repeat the calculation of det αi and det βi Nt times. Therefore, the complexity of its
29

()
3 × N . On the other hand, the complexity of naive
algorithm is O
(N
N
N
×
4
×
3)
x
y
z
t
(
)
method is O (Nx Ny Nz Nt × 4 × 3)3 . If we set Nx = Ny = Nz = Nt = 10, it is reduced
from O(1015 ) to O(1013 ). Then, the running time is decreased to 1/100, roughly.
When we use the order function O(·), we usually neglect the constants. In following
sections, we only consider the dependence of Nx , Ny , Nz and Nt because the physics ofsystem depends on them.
Although the reduction formula also seems eﬃcient, calculation cost is still large.
Especially,( when we3 )calculate with large lattice, the time complexity becomes to diverge
with O (Nx Ny Nz ) . Therefore, we need to reduce it.

6.2

Hopping parameter expansion

To reduce the time complexity of the determinant calculation, we can use the hopping
parameter expansion (HPE). In this method, we firstly rewrite the determinant as follows:[]det ∆ = exp Tr log ∆
(88)
[
(
)]
= exp Tr log 1 − κQ .
(89)
Then we use the fact that the Wilson fermions form as ∆ = 1 − κQ.
Here, if hopping parameter κ is small enough, we can expand the logarithm as
[
]
∞
∑
1
det ∆(µ) = exp Tr(κQ(µ))k .
(90)
k
k=1

This is the hopping parameter expansion. Using this equation,
we can calculate
det ∆(µ)()
by the polynomial of Q(µ). Then, the complexity is O (Nx Ny Nz Nt )2 . This is because
the trace of Qk is calculated by multiplication of Q and a vector. Such computation2 ) operations. We can reduce the time complexity using this method.
costs us O(Nrank
However, the calculation cost is still large. It is because we have to calculate det ∆(µ)
for varying µ. In order to reduce the cost, we can use the following method.

6.3

Winding number expansion

Winding number expansion (WNE)[35, 40] is an improved method of HPE.
In HPE, we rewrite det ∆ to exp [Tr log ∆] and expand log ∆. Then, using a methodof statistical mechanics, we can rewrite the trace as∑⟨⟩
Tr log ∆(µ) =Ψ log ∆(µ) Ψ ,(91)
Ψ

where |Ψ⟩ is a state in the system. The right hand side of this equation means asummation of paths of which the initial state and the final state are the same (closed
loops). Thus, we can rewrite the trace to a summation of closed loops.

30Here, there are many types of size or shape of loops. We therefore have to distinguishthem to calculate in practice. Now, let us focus on the chemical potential dependence.
From the definition of ∆(µ), we notice that the coeﬃcient exp(µaNt ) = exp(µ/T )appears when a fermion hop Nt times along the time direction. We hence can write
det ∆(µ) formally,
[ ∞
]
)k
(
∑
det ∆(µ) = exp
Wk eµ/T
.
(92)
k=−∞

We can calculate Wk numerically. The algorithm is discussed in Appxs. B and C. (Seealso Refs. [41] and [42].) Using this equation, we can see the µ dependence of det ∆clearly.
Note that Eq. (92) can be used at any chemical potential. That is, we can calculate Z(µ) from det ∆(µ) via WNE without canonical ensembles. Such a calculation is
discussed in Appx. D.
Since Wk does not depend on µ, if ONCE
( we measure Wk , we
) can calculate det ∆with arbitrary µ. The time complexity is O (Nx Ny Nz )2 × (Nt )2 in this method. (SeeAppxs. B and C.) The complexity of spatial part of this method is smaller than the oneof reduction formula. Therefore, we can calculate ZC (n) in large lattice using WNE. In
this thesis, we adopt this method when we calculate ZC (n).

7Problems of canonical approach

It is known that canonical approach has several numerical problems. In this section, we
will see the problems using a raw data of canonical partition function as an example.
The following figure shows a raw data of canonical partition function calculated on(1)
one configuration ZC (n). We then generate a configuration at T /Tc = 0.93(5) andmπ /mρ = 0.80, and set the sample size of Fourier transform as NF T = 1024. We findthat ZC (n) is a monotonic decrease function when n is small and it becomes smaller
rapidly with increasing n. However,• ZC (n) is saturated at n = 50, and• it recovers near n = 1023 = NF T − 1.
As we know, ZC (n) describes an existence probability of a system which n particles are
in. Therefore, above figure means that, when µ = 0, n = 1 system and n = 1023 system
exists with the same probability!
It is manifestly ill data. In following sections, we discuss why such data were generated.

7.1

Property of discrete Fourier transform

First, we remember a property of discrete Fourier transform (DFT).

31

Figure 16: Raw data of canonical partition function with one configuration.
DFT is a transformation which is defined as the following equation:F̃ (k) =

1
NF T

N∑
F T −1f (j)e−iϕk j ,

ϕk =

j=0

2πk
.
NF T

(93)

Now, we only have f (j) (j ∈ Z, 0 ≤ j < NF T ) which are sampled on NF T points. In
DFT, we construct the finite series F̃ (k) using the finite series f (j).
Then, we note that,exp (−i ϕk j) = exp (−i [ϕk + 2πn] j) = exp (−i [ϕk+NF T n ] j) ,(94)

with any n ∈ Z. That is why we can not distinguish F̃ (k) from F̃ (k + NF T n). Thisphenomenon is called aliasing.
Aliasing comes from the periodic boundary condition of DFT, and it makes us equate
low negative frequency with high positive frequency. We can see this considering n = 1
and k ≤ 0 case for example. Concretely,F̃ (0) = F̃ (NF T )(95)F̃ (−1) = F̃ (NF T − 1)

(96)

F̃ (−2) = F̃ (NF T − 2) .

(97)

Thus, Fig. 16 shows that n = 1 system and “n = −1” system exists with the sameprobability when µ = 0. This means that the symmetry of the exchange of particles foranti–particles (C–symmetry) is satisfied in the QCD system. Then, the highest baryon
number is nmax = NF T /2.
32

Such a property is seen in the work of H. Nyquist[43] in 1928 and proved by C.E.
Shannon[44] in 1949. nmax = NF T /2 is called Nyquist frequency or the folding frequency.
We may see this property from Fig. 16. We noticed that ZC (n) in Fig. 16 is symmetricaround n = NF T /2 = 512. Figure 17 shows it clearly. Note that this data is not

Figure 17: Raw data of canonical partition function with one configuration. Blue line
denotes half of data.
symmetric completely because there is the numerical errors which we will discuss below
in detail.
Therefore, although we calculate ZC (n) using DFT with NF T , we can use only halfof data.

7.2∗

An instability of DFT

We do not understand yet why ZC (n) is saturated at n = 50. Although it had notbeen investigated for a long time, the author studied it from a numerical standpoint
recently[45, 46].
First, the author considered that this saturation is caused by the numerical errors,
and if we understand how the errors aﬀect to our calculation, we can reduce this saturation. He notes that we can classify the numerical errors into four types of errors[45](see also Refs.[47] and [48]).
1. Rounding error: for a numerical calculation, numerical values are always rounded
when they are assigned to each variable because the variables have only finite
significant digits. Therefore, we can not represent a mathematical value completely
by using the numerical variables; the diﬀerence of mathematical and numericalvalues is called rounding error.
33

For example, this error arises in the following type of calculation,
1= 3.333333 · · · × 10−1 (Exact mathematical value)3
assignment−−−−−−−→ 3.333333 × 10−1 (Numerical value; 7 digits.) ,(98)

In this case, the exact mathematical value has infinite significant digits, but thenumerical value has only seven significant digits. The numerical result is less thanthe exact mathematical one as 3.333 · · · × 10−8 . To perform rounding exactly,
computers usually have a guard digit[48].
2. Truncation error: when we estimate an infinite series numerically, we have to
truncate it because we can use only finite computation resource. Then, there isa diﬀerence between the exact mathematical result and its truncated one. Thisdiﬀerence is called truncation error. We can estimate the eﬀect of this error by
varying a truncation point of the series.
3. Cancellation of significant digits: when we calculate a subtraction between two
nearly equal values, a lot of higher significant digits are canceled. Therefore, theresult has only a few significant digits because the number of significant digits islimited in the numerical calculation. This phenomenon is called the cancellationof significant digits.
It occurs in the following type of calculation:
1.234567 − 1.234566 (7 digits)= 0.000001 (1 digit.) .

(99)

In this case, variables have seven significant digits at first. However, six significant
digits are lost when we calculate the subtraction. If we want to reduce the eﬀect of
this cancellation, we should increase the number of significant digits. For instance,
we can consider the following calculation instead of above one:
1.234567444444444444444 − 1.234566111111111111111 (22 digits)= 0.000001333333333333333 (16 digits.)(100)

Then, although six significant digits are similarly lost in this calculation, sixteen
digits still remain in the final result.
4. Loss of trailing digits: when we calculate an addition or subtraction between ahuge number and small one, many lower significant digits of the small number are
cut oﬀ. This is because the variables hold only finite significant digits. Then, thenumerical result is underestimated. This phenomenon is called the loss of trailing
digits.

34

Specifically, it arises the following type of calculation:
3 × 1010 + 2 × 100= 3.000000002 × 1010 (Exact mathematical result)assignment−−−−−−−→ 3.000000 × 1010 (Numerical result; 7 digits.) ,(101)

In this case, variables have seven significant digits, and the numerical result isunderestimated from the mathematical one by 2 × 100 .
As the cancellation of significant digits, we should increase the number of significant
digits in order to reduce this loss.
3 × 1010 + 2 × 100= 3.000000002 × 1010 (Exact mathematical result)assignment−−−−−−−→ 3.0000002000000 × 1010 (Numerical result; 16 digits.) ,(102)

In this calculation, by increasing significant digits, we reduce the error between
the mathematical result and numerical one.
Next, the author considered as follows.
• Usually, the rounding error is not transmitted to higher significant digits sincecomputers have a guard digit[48]. Sometimes we may see that this error transmits
to higher significant digits, but actually, only the cancellation of significant digits
or loss of trailing digits occurs at the same time. In this paper, we consider thetransmission of rounding error as an eﬀect of the cancellation or loss.
• The truncation error may occur in our program. When we use HPE or WNE inorder to calculate ∆(µ), we have to truncate the infinite series of these expansion.
Although the truncation error arises in our calculation, we can control it andestimate its eﬀect by changing a truncation point of this series.
Thus, we conclude that the instability which we do not control is caused by the cancellation of significant digits or the loss of trailing digits, or both.
Actually, he monitored the behavior of all variables in a DFT program, and foundthat the cancellation of significant digits occurs and cause the instability of DFT. Thiscancellation arises almost at the last addition of the summation in Eq. (93); specifically,[(N −2 (
)
)
∑
1
2πk
−i(2πk/N )ne
ZC (n) =
Z
N
N
k=0
]
(
)
2π(N − 1) −i(2π(N −1)/N )n
+Z
e
,
(103)
N
↑ This addition.
35

Figure 18: Number of canceled significant digits as a function of baryon number. The
red and green points are data of low temperature and high temperature.
Figure 18 shows how many digits are lost at the summation of Eq. (103). (Note thatthe author then uses the multi precision calculation which is described in section below.)
The algorithm is discussed in detail in Appx. E.
From this figure, we can see that this cancellation becomes serious with increasingn. For example, about 60 digits are canceled at n = 100, 130 digits are canceled at
n = 200. Usually, one use the double precision variables in order to calculate thephysical quantities, but it has only 16 significant digits. Thus, we can not calculate
ZC (n) well using the usual (double precision) variables.

7.3∗

Reduce these problems

In order to reduce these problems, the author suggests to calculate ZC (n) with large
enough data space, specifically,
1. NF T should be set large enough value for one’s purpose.
For instance, the author aims to study the finite density phase transition of QCD,
i.e. the confinement–deconfinement phase transition. We will use a small lattice inthe discussion below, whose size is 83 × 4. One considers that QCD critical point
may be at roughly Tc ≈ 160 MeV and µc = µcB /3 ≈ 200–400 MeV[12, 13, 14].
Therefore, at T = Tc , lattice spacing a is about 0.3 fm because Nt a = 1/T in finite
temperature lattice QCD. Then, the lattice size is (2.4 fm)3 × (1.0 × 10−24 s).
If critical chemical potential µc is 200 MeV, it corresponds to length l = 1.0
fm. We then may consider that l means the distance between quarks. Thus, by
putting about 10–20 quarks in the system, we may simulate near the critical point.
Accordingly, we need to set NF T /2 > 20.
36

2. The number of significant digits should be set too.
The magnitude of the cancellation of significant digits is related with the net
quark number n. If NF T is large, the cancellation is serious, and then we have toincrease the number of significant digits. (We are going to see how to increase it
below.)
The canceled significant digits are roughly proportional to n. Therefore, one has toestimate or measure how many digits are canceled and set the number of significant
digits to large enough value.
The author studied how to enhance the accuracy of ZC (n) using the multi precisioncalculation. This calculation extends the number of significant digits. We can make itby binding several variables. For instance, a real number x represents as follows in thebase–b numerical system,()
x = xsign · x1 b−1 + x2 b−2 + · · · + xf b−f · bx0 ,
(104)
where
• xsign is the signature of x (xsign ∈ {+1, −1});• x0 is the exponent of x (x0 ∈ Z);• x1 , x2 , · · · , xf are the mantissas (or significands, x1 ∈ {1, 2, · · · , b − 1}, xi ∈{0, 1, 2, · · · , b − 1} (i = 2, 3, · · · , f ));• f is the number of significant digits; and• b is the base (or radix).
Usually, the computers are designed with b = 2 and fixed f . Then, xi are just the bits,
x0 is a set of several bits. However, if we program this framework as a software (xi areconstructed by integer variables with any f , also b is a variable), we can change b andf , and then we can use a original variable with wide enough range and high accuracy.
In this thesis, we use the libraries of multi precision calculation, MPFR (C++) and
FMlib (FORTRAN) for our simulation.
We can reduce the cancellation error by increasing the number of significant digits
with multi precision. The author calculated ZC (n) changing the precision of variables,
and compared them as Fig. 19 shows. We find that, by increasing the number of significant digits f , we can calculate ZC (n) up to large n. This is because ZC (n) is justa remainder of the calculation of Eq. (103). That is, the number of canceled significant
digits is directly related with a magnitude of ZC (n).
He notes that, when we change the lattice volume V or temperature T fixing thenumber of significant digits, a saturation point (kink) of ZC (n) moves.
ZC (n) describes the existence probability of a system which n particles are in. If
only T becomes large, excited states become easy to appear in the system, then ZC (n)decreases slowly with increasing n. We can see this actually in Fig. 19.
37

Figure 19: Normalized canonical partition function as a function of baryon number. The
left and right panels show data of low temperature and high temperature, respectively.
The red, green, blue and cyan points are the data which is calculated 16, 32, 48 and 64
significant digits, respectively.
In the other situation, if only V becomes large, also ZC (n) decreases slowly too since
the net quark number density n/V becomes small. If volume is changed as V → V ′ = 2V ,ZC (n) corresponds to ZC′ (n′ = 2n) because they represent the same density system.
Thus, we can replace ZC (1) → ZC′ (2), ZC (2) → ZC′ (4), · · · , therefore ZC (n) decreases
slowly.

8

Study of finite density phase transition with canonicalapproach

As seen above, by introducing these methods, we can avoid the sign problem and calculate ZC (n) with high accuracy. In this thesis, we calculate several thermodynamic
observables using the canonical approach with the multi precision calculation. We will
see the results below in detail.

8.1

Thermodynamical observables and partition function

From grand partition function Z(µ), we can calculate thermodynamical observables asfollows,p(µB , T )1=
log Z(µB , T )
T4
V T3
Nt 3
log Z(µB , T ) ,=Nx Ny Nz(105)

where p is the pressure, V = Nx a × Ny a × Nz a is the spatial volume of the system
and µB = 3µ is the baryon chemical potential. We then use a relation in the finite
temperature lattice QCD, Nt a = 1/T , with lattice spacing a.
38

Taking a deviation of the pressure, we can get a useful relation,∆p(µB , T )p(µB , T ) p(0, T )
=
−
4
4
T
T4
(T
)
3
Nt
Z(µB , T )=log.
Nx Ny Nz
Z(0, T )(106)

It is useful because the partition function is normalized by Z(0, T ) as the expectation
values of observables.
Other observables, baryon number density nB and baryon number susceptibility χ
are given by the diﬀerential of pressure,nB (µB , T )∂p(µB , T )=,3T∂(µB /T ) T 4
χ(µB , T )∂2p(µB , T )
=
.
2
2
T
∂(µB /T )T4

Note that, the partition function has the following periodicity as Eq. (20),
(
)
(µ)
µ
2πk
Z
=Z
+i
,
T
T
3

(107)
(108)

(109)

which is caused by RW symmetry. Using this equation, we can rewrite the partition
function as[ 2()](µ) 1 ∑
µ
2πk
=
Z
+i
.
(110)
Z
T
3
T
3
k=0

Then, the canonical partition function ZC (n) becomes()]∫ 2π ( ) [ ∑
2
1
1
µI
iµI
2πk
ZC (n, T ) =dZ+iei(µI /T )n
2π 0
T
3
T
3
k=0
(
)]
∫ 2π ( ) [1µI
1 + e2πi/3 + e4πi/3
iµI
=
d
Z
ei(µI /T )n .
2π 0
T
3
T

(111)

Therefore, ZC has the following property,
ZC (n ̸= 3k) = 0 ,(112)

for any k ∈ Z. This property is called the triality.
Thus, the grand partition function is written as
Z(µB ) =

∞
∑

(
)B
,
ZC (B) eµB /T

B=−∞

39

(113)

where B = n/3 is the baryon number. Because of this, we use the baryon number B
and baryon chemical potential µB instead of the quark number n and quark chemical
potential µ below.
In the following section, we compare our results of the canonical approach with theresults of MPR in order to discuss the validity range of canonical approach. We use theresults of MPR in Ref. [24], and we adopt the same numerical setup as Ref. [24] when
we calculate our observables of canonical approach.
In Ref.[24], the authors also discuss consistency between MPR and Taylor expansion
method with Wilson–clover fermions and concluded that both methods produced consistent results in a small chemical potential region where errors of both methods could beunder control. Therefore, our work enables us to check consistency among our canonicalapproach, MPR and Taylor expansion method.

8.2

Numerical setup

We adopted the clover improved Wilson fermion action with Nf = 2 and CSW =(1 − 0.8412/β)−3/4 evaluated by a one–loop perturbation theory[17], and Iwasaki gauge
action[16]. All simulations were performed on a Nx × Ny × Nz × Nt = 8 × 8 × 8 × 4
lattice. Number of flavor is set to Nf = 2. We considered values of β = 2.00, 1.95, 1.90,
1.85, 1.80 and 1.70; they correspond to T /Tc = 1.35(7), 1.20(6), 1.08(5), 0.99(5), 0.93(5)and 0.84(4), respectively. The values of the hopping parameter κ were determined foreach value of β by following the line of constant physics for the case of mπ /mρ = 0.80,
as in Ref. [28]. In Ref.[28], the lined of constant T /Tc is determined by the relation,
Tc /mρ = A(1 + B(mπ /mρ )2 )/(1 + C(mπ /mρ )2 ) with A = 0.2253(71), B = −0.933(17)and C = −0.820(39), which obtained in Ref. [49] to evaluate Tc /mρ for each mπ /mρ .
We generated gauge configurations at µ = 0 using the HMC method. We then set
the step size dτ and number of steps Nτ for HMD part to dτ = 0.2 and Nτ = 50 sothat the simulation time was Nτ × dτ = 1. We neglected the first 2000 trajectories
for termalization, after we extracted 400 configurations every 200 trajectory for eachparameter set.
Number of Fourier series NF T is set to NF T = 1024 for test in this thesis. When we
calculate the grand partition function Z(µ) from the canonical partition function ZC (n),we have to truncate the series of fugacity expansion. The truncation point Bmax is listed
in Tab. 2, which is discussed below in detail. Number of significant digits f is set to
f = 400 in this thesis.

8.3

Error estimation

When we calculate the fugacity expansion of the grand partition function Z(µ), we musttruncate the series as
Z(µB ) =()B.
ZC (B) eµB /T

B
max∑
B=−Bmax

40

(114)

β
2.00
1.95
1.90
1.85
1.80
1.70

CSW
1.50579
1.527171.55044
1.57589
1.60380
1.66885κ
0.136931
0.137716
0.138817
0.140070
0.141139
0.142871

T /Tc
1.35(7)1.20(6)1.08(5)0.99(5)0.93(5)0.84(4)mπ /mρ
0.80
0.80
0.80
0.80
0.80
0.80

# of Conf.
400
400
400
400
400
400

Bmax
36
34
28
34
12
30

Table 2: Parameters of each simulation.
Thus, the truncation error arises, and then we have to estimate how large the error is.
In order to estimate it, we vary the number of series Bmax and check the convergenceof results. If the results are stable when Bmax is varied, we can trust them because theyshould not depend on Bmax physically.
For a convergence check, we evaluate the relative error R as follows[46, 50],
R(µB ) :=⟨O(µB )⟩Bmax − ⟨O(µB )⟩Bmax −1,⟨O(µB )⟩Bmax(115)

where ⟨O(µB )⟩Bmax is the expectation value with the truncation at Bmax . In this thesis,
we require that the results have at least two significant digits. Therefore, we trust theresults only forR(µB ) < 10−3 .
(116)

8.4

∗

Numerical results

Using the error estimation method, we study the chemical potential dependence of thethermodynamic observables and estimate the validity range of our canonical approach.
8.4.1

Canonical partition function and analysis of truncation point

Before we introduce our results of thermodynamic observables, we show the results ofcanonical partition function. Figure 20 shows the baryon number dependence of thecanonical partition function. From this figure, we can see that ZC (B) becomes smaller
rapidly with increasing B at both of temperatures below and above Tc . Also, we seethat ZC (B) depends on temperature; it is raised with increasing temperature. We canconsider that they are because ZC (B) represents the existence probability of the system
which B particles are in*3 . That is, B = 1 system is more likely than B = 2 or anothersystem when µ = 0, and B = 2 system at high temperature is easier to arise than such
systems at low temperature.
In Fig. 20, the canonical partition function is cut at Bmax . We define the value of Bmax
using the following way. Because of numerical and statistical errors, sometimes ZC (B)becomes negative, although it has to satisfy ZC (B) > 0. Then, if the error of ZC (B) is*3
In detail, ZC (B)(eµ/T )B /Z(µ) roles the weight function as we saw in Sec. 5, and then it represents
the existence probability of the system. When µ = 0, ZC (B)/Z(0) represents the probability. We now
neglect the factor 1/Z(0) in the discussion of this section because the factor is constant for B.

41

Figure 20: Baryon number dependence of normalized canonical partition function.
larger than the absolute value of ZC (B), we can use the positive value −ZC (B) insteadof the negative value ZC (B). However, if ZC (B) is negative and δ [ZC (B)] < |ZC (B)|,
we have to consider that the calculation does not work well at this parameter. Therefore,
when we find such ill data, we reject this and later data. Then, we set the truncation
point Bmax from the previous data.
8.4.2

Pressure, baryon number density and baryon number susceptibility

Let us introduce our results of thermodynamic observables.
First, we examine the pressure. Fig. 21 shows the chemical potential dependence ofpressure. From this figure, it is found that the results of pressure at high temperature
T > Tc do not suﬀer from large errors up to µB /T = 5, approximately. The results at
low temperature T < Tc are reliable up to µB /T = 3.5–4. Conversely, they are reliable
only up to µB /T ≈ 3. This may be because we generated configurations at µ0 = 0
and they suﬀered from fluctuations caused by the (crossover) phase transition at zerodensity. We may obtain clearer signals if we generate configurations at pure imaginary
chemical potentials because we can leave from the crossover phase transition line.
Fig. 22 shows the comparison of pressure calculated by the canonical approach and
MPR method. We can see that the results of canonical approach are consistent with theresults of MPR method.
Next, we consider the expectation value of the baryon number density. Fig. 23 shows
the chemical potential dependence of baryon number density. This figure demonstrates
42

Figure 21: Chemical potential dependence of pressure. The red, green, blue, cyan,
magenta, and brown points are the results at T /Tc = 1.35(7), 1.20(6), 1.08(5), 0.99(5),0.93(5) and 0.84(4), respectively. The upper bound of the baryon chemical potential isdetermined by Eq. (116). This figure was already published in Ref. [45]that, for temperatures above and below Tc , the results are reliable up to µB /T = 4 andµB /T = 3–3.5, respectively. However, near Tc , the reliable chemical potential range arelimited up to µB /T = 2.4. This may be for the same reason described in the pressure
analysis.
Figure 24 represents the comparison of baryon number density between the methods.
It demonstrates good agreement between the results of the canonical approach and thoseof MPR method in this case.
Moreover, the gradient of the baryon number density nB /T 3 as a function of baryon
chemical potential µB becomes small as the temperature decreases. In a zero temperaturecase, nB is expected to be zero up to µB = mB where mB is the lightest baryon massof a system, and becomes a finite value at this point. This phenomenon is called thesilver blaze phenomenon[51]. The data at T /Tc = 0.84 in Fig. 23 does in fact show sucha feature.
Finally, we investigate the baryon susceptibility. Figure 25 shows the results of thesusceptibility. From this figure, we find that the results are reliable up to µB /T = 3.5
at T > Tc , whereas they are reliable up to Tc = 2.4–2.9 at T < Tc .
Figure. 26 represents the comparison between the methods in the susceptibility case.
We see that the susceptibility results of canonical approach are in very good agreementwith those of the MPR method.
The baryon susceptibility as a function of µB /T does not show a clear peak. Thus,
43Figure 22: Comparison of pressure calculated by the canonical approach and the MPR
method. The colors of the data points are the same as in Fig. 21 with several additional
colors. The data points plotted in the additional colors of dark red, dark green, dark
blue, dark cyan, dark magenta, and dark brown points are the results at T /Tc = 1.35(7),1.20(6), 1.08(5), 0.99(5), 0.93(5) and 0.84(4), respectively, as calculated by the MPR
method. This figure was already published in Ref. [45]Figure 23: Chemical potential dependence of baryon number density. The colors ofthe data points are the same as in Fig. 21. The upper bound of the baryon chemical
potential is determined by Eq. (116). This figure was already published in Ref. [45]44Figure 24: Comparison of the baryon number densities calculated by the canonicalapproach and the MPR method. The colors of the data points are the same as inFig. 22. This figure was already published in Ref. [45]Figure 25: Chemical potential dependence of baryon susceptibility. The colors of thedata points are the same as in Fig. 21. The upper bound of the baryon chemical potential
is determined by Eq. (116). This figure was already published in Ref. [45]45Figure 26: Comparison of baryon susceptibilities calculated by the canonical approachand the MPR method. The colors of the data points are the same as in Fig. 22. This
figure was already published in Ref. [45]signals of the finite density transition between the confinement–deconfinement phases is
not observed yet in this case.

46

9Conclusion and discussion

In this thesis, we have shown that the canonical approach is consistent with MPR
method. Moreover, the canonical approach provides reliable results beyond µB /T = 3
for almost all observables. This is very encouraging for the first–principles calculation offinite density QCD because other methods, such as MPR, Taylor expansion and imaginary chemical potential method[52], yield reliable information in practical situations
only up to µB /T = 3. Here, the multi precision calculation significantly contributes tothis conclusion.
We have investigated the QCD phase structure via the pressure, baryon number density and baryon number susceptibility. We calculated these thermodynamic observablesin this thesis. Our results have not shown the QCD phase transition, yet. However,
it is important that we reduced the instability of canonical approach using the multi
precision calculation and calculated the observables with higher accuracy than previous
studies.
In order to get more reliable signals of thermodynamic quantities in a large baryon
chemical potential region, we need to calculate the canonical partition functions moreaccurately at large baryon numbers. As shown in Fig. 19, our canonical partition functions currently have unphysical phases at several baryon numbers although these shouldhave a real positive value in principle. After the author works, ones have studied whysuch phases appear in the canonical partition function ZC (n)[53, 54]. However, we havenot understood it well.
In our opinion, this may be originated from the overlap problem[18]. In this work, we
calculate the grand canonical partition functions at pure imaginary chemical potential
with gauge configurations generated at zero chemical potential through the simplest
reweighting method. However, practically speaking, we have to calculate them withgauge configurations generated at each pure imaginary chemical potentials to realize theappropriate importance sampling. This is a weak point of our work and we need to
improve it in future.
After the author’s study was published, the improvement of canonical approach havebeen applied to the other works. For example, the author and his collaborators calculated
the higher cumulants with our improvement of canonical approach[55]. However, theerror bars are still large in this study because of the weak point that we discussed above.
On the other hand, D.L. Boyda, et al. (the lattice QCD group of Vladivostok) used
our improvement and studied how to calculate the canonical partition function ZC (n)with more high accuracy[56, 57]. Then, they used fitting procedure for the number
density at imaginary chemical potential and constructed ZC (n) from it with more highaccuracy.
In addition, they calculated the higher cumulants with their procedure based on
our improvement of canonical approach[58]. Then, they showed that their results wasconsistent to experimental results.

47

The author notes that previous studies of sign problem are the mere resummation
techniques of a series or gauge configurations. In the reweighting technique, we extract
the phase of probability and resum up gauge configurations with new probability. In the
Taylor expansion method and canonical approach, we expand the observables using the
Taylor, fugacity and winding number expansion. In addition, recently ones study the
Lefschetz thimble to reduce the sign problem. This is just a resummation technique ofgauge configurations.
The author considers that such methods may include the calculation which causes thecancellation of significant digits. This is because, in the canonical approach, we can not
avoid the cancellation even though we use and resum the discrete Fourier transformation
or fast Fourier transformation. Also he consider that the multi precision calculation isalways eﬀectual for such methods.
The author notices that we can calculate ZC (n) with high enough accuracy using
the reduction formula even though we do not use the multi precision calculation*4*5 .
This means that, in the reduction formula, small values result from the multiplication,
not from the subtraction between two large values. From the algorithm of reduction
formula, we may study how to reduce the cancellation of significant digits.
The canonical approach has been investigated in the previous studies [29, 32, 33,
34, 35, 36, 59, 60, 61, 62, 63]. We can also find through our work that the canonicalapproach is a useful and promising method. However, our method may be improved
further to obtain results under more realistic conditions, i.e. a lighter quark mass, larger
volume, a finer lattice spacing and higher density. Although the hopping parameter
expansion yielded very interesting results in this thesis, we want to calculate the fermion
determinant without the approximation which is from the truncation, We have learned
from this study that the key point is how to calculate the determinant at imaginary
chemical potential values in order to compute the Fourier transformation in Eq. (67)with high accuracy. This requires more computational resources than what has beenreported here but is within the scope of the next–generation high–performance era.

*4Then, we only need to compute the mantissas and exponent part respectively. This is because verysmall eigenvalues appear (λ < 10−308 ) in this calculation. We do not need to extend the number ofsignificant digits.
*5But, as seen above, the calculation cost of reduction formula is much larger than the cost of winding
number expansion and Fourier transformation.

48

Acknowledgment
I am deeply indebted to my supervisors, Prof. Tohru Eguchi and Prof. Hidekazu Tanaka
for insightful comments and beneficial support. In particular, my laptop which I am
using right now was given by Prof. Eguchi. Also I would like to thank my sub–chief
examiners, Prof. Yu Nakayama and Prof. Jiro Murata for useful advice and comments.
I am greatly thankful to my collaborators, especially Prof. Atsushi Nakamura and
Dr. Ryutaro Fukuda, for productive discussion, significant support and sincere encouragement. In particular, my tablet PC which I had used for study had been given by
Prof. Nakamura. I am grateful to Mr. Ryo Sugawara for linguistic support.
I would like to thank everyone who has supported me, especially those who gave mebeer, sake, wine, comics, PCs and money.
My calculations were performed using the supercomputers, SX–9 and SX–ACE atReserach Center for Nuclear Physics (Osaka Univ.), PACS–IX and HA–PACS at Centerfor Computational Sciences (Univ. of Tsukuba), SR16000 at Yukawa Institute for Theoretical Physics (Kyoto Univ.) and FX10 at Information Technology Center (Univ. ofTokyo).
My study was partially supported by the Large Scale Simulation Program in KEK[No.14/15-19 (FY2014)], Interdisciplinary Computational Science Program in CCS [16a44],
joint research program in JHPCN [EX1631], HPCI System Research Project [hp170197]
and KAKENHI [26610072, 24340054, 22540265, 15H03663 and 26610072].

49

Appendix
A

Grassmann integral of fermion action

Fermion fields are the Grassmann numbers which satisfy the anticommutation relation,
Ψi Ψj = −Ψj Ψi ,(117)

Ψ̄i Ψ̄j = −Ψ̄j Ψ̄i .

(118)

The Grasmann integral is defined by the following relations,∫∫dΨi = dΨ̄i = 0 ,∫∫dΨi Ψi = dΨ̄i Ψ̄i = 1 .

(119)
(120)

Also, we get the following relations from Eq. (117),
Ψi Ψi = Ψ̄i Ψ̄i = 0 .

(121)

From this equation, we find that the polynomial of Grassmann numbers includes theterms of Ψi and Ψ̄i finitely.
We perform the Grassmann integral theoretically using these equation. For simplification, we consider two flavor fermions Ψi (i = 1, 2) and their chiral partners Ψ̄i(i = 1, 2). The integral of exp(−Ψ̄AΨ) (A is a matrix) is calculated as∫∑2dΨ1 dΨ̄1 dΨ2 dΨ̄2 e− i,j=1 Ψ̄i Aij Ψj∫= dΨ1 dΨ̄1 dΨ2 dΨ̄2 e−Ψ̄1 A11 Ψ1 e−Ψ̄1 A12 Ψ2 e−Ψ̄2 A21 Ψ1 e−Ψ̄2 A22 Ψ2∫()()= dΨ1 dΨ̄1 dΨ2 dΨ̄2 1 − Ψ̄1 A11 Ψ1 1 − Ψ̄1 A12 Ψ2()()
× 1 − Ψ̄2 A21 Ψ1 1 − Ψ̄2 A22 Ψ2= (−A11 )(−A22 ) + (−A21 )(+A12 ) = det A .

(122)

Then, we use the Taylor expansion of the exponential in the third line. Note that thehigher powers (e.g. (Ψ̄1 A11 Ψ1 )2 ) do not appear because of Eq. (121). We also use
Eq. (119) and Eq. (120) in the fourth line.
We can generalize this equation to more general formulation. T. Matthews and A.
Salam give us the following formula[64],∫DΨDΨ̄ e−Ψ̄∆Ψ = det ∆ ,
(123)
∫
(
)
(124)
DΨDΨ̄ Ψi Ψ̄j e−Ψ̄∆Ψ = − ∆−1 ij det ∆ ,where ∆ is the fermion matrix. This is called Matthews–Salam formula.
50

B

∗

Algorithm of winding number expansion

We can calculate the coeﬃcient Wk of Eq. (92) via HPE[42, 41]. This is because TrQm
is given as follows,∑
TrQm =⟨Ψa,α (x) |Qm | Ψa,α (x)⟩ ,
(125)
x,a,α

where a and α denote the Dirac and color index, respectively. Here, we dissolve Q into(+)
the spatial component Qi (i = 1, 2, 3), forward temporal component Q4 and backward(−)
temporal component Q4 ,
( 3
)m
∑
m
a,α
+µa (+)−µa (−)
Q |Ψ (x)⟩ =Qi + eQ4 + eQ4|Ψa,α (x)⟩ .
(126)
i=1

Note that Qi includes the clover term. Then, we can classify this with chemical potential
dependence. That is,+m∑(

Qm |Ψa,α (x)⟩ =

k=−m
+m
∑

=:

e+µa

)k

Q(k) (m) |Ψa,α (x)⟩( +µa )k (k)⃗ (m; x, a, α) ,
e
X

(127)

k=−m

where Q(k) (m) is the product of Qi and Q±
4.
m
Thus, we can write TrQ as
⟨
⟩
+m
∑
∑
( +µa )k (k)m⃗⃗ (m; x, a, α) ,
TrQ =
Ψ(x,
a, α),
e
X
x,a,α

(128)

k=−m

⃗ := |Ψ⟩ and ⟨·, ·⟩ is the inner product.
where ΨHere, note that k denotes the time slice because it is the power of eµa . Then, we
can neglect the variable t of the vectors,
therefore we can) reduce( the time complexity()⃗
of the multiplication QΨ from O (Nx Ny Nz Nt × 3 × 4)2 to O (Nx Ny Nz × 3 × 4)2 .
In⟨ the practical⟩ computing, we do not perform the summation of m. We thus save
µa
k
(k)
⃗
⃗
Ψ, (e ) X=: Z (k) as an array for k. After we perform the summation of m, we
extract the coeﬃcient of (eµ/T )k = (eµa )kNt . The coeﬃcient is just Wk .
Because the algorithm is intricate, we show it as a pseudo code on this page. Notethat the outside loops of (x, y, z, t, c, µ) for the trace are performed by the noise method
which is discussed below.

C

Noise method

The trace of matrix is performed by the loops of summation. Naively, we perform theloops Nx Ny Nz Nt ×3×4 times. We can reduce this calculation using the noise method[65].
51

Algorithm 1 Algorithm of WNE via HPE
Initialize: x, y, z, t, c, µ, k, m ∈ Z
Initialize: X (k) (x, y, z, c, µ), Y (k) (x, y, z, c, µ) ∈ C
Initialize: Z (k) ∈ C
Initialize: S(x, y, z, c, µ), T (x, y, z, c, µ) ∈ C
Initialize: W (k) ∈ C
for (x, y, z, t, c, µ) = (1, 1, 1, 1, 1, 1) to (Nx , Ny , Nz , Nt , 3, 4) dofor k = −kmax to +kmax do
Z (k) = 0
⃗ (k) = ⃗0X
end forX (0) (x, y, z, c, µ) = 1
for m = 1 to mmax dokrange = m
for k = −krange to +krange do⃗ = κ ∗ (Q1 (t) + Q2 (t) + Q3 (t)) ∗ X
⃗ (k)
S
(+)
(−)
⃗ (k−1) + Q (t) ∗ X⃗ (k+1) )T⃗ = κ ∗ (Q4 (t) ∗ X
4(k)⃗⃗ + T⃗Y=S
end forfor k = −krange to +krange do
Z (k) = Z (k) + Y (k) (x, y, z, c, µ)/m
⃗ (k) = Y⃗ (k)X
end forend forfor k = −(krange − 1) to +(krange − 1) do
if mod(k, Nt ) ̸= 0 thenW (k/Nt ) = −Z (k)end ifend forend for52

In this method, we start from Eq. (13), that is,∑⟨ϕ |A| ϕ⟩ ,
TrA =

(129)

ϕ

where A is a matrix with the rank of Nrank . Then, |ϕ⟩ is a vector. Here, we use the⃗ := |ϕ⟩ as above section, we then write as
notation ϕ
⟩
∑⟨
⃗ Aϕ
⃗ ,
TrA =
ϕ,
(130)
ϕ

where ⟨·, ·⟩ is the inner product.
We can calculate it choosing ϕ as the vector of noises[65],⟩1 ∑ ⟨ (i)TrA = lim
⃗η , A⃗η (i) ,N →∞ N
N

(131)

i=1

where ⃗η (i) is the i–th noise vector. The noise vector satisfies the following property,
1 ⟨ (i) (j) ⟩⃗η , ⃗η= δij .
(132)
Nrank
⃗
This method corresponds to the Monte Carlo method for each component of ϕ.
J. Foley, et al. showed that we can approximate the trace by 16 noise vectors withhigh enough accuracy[65]. In this thesis, for safety, we use 64 noise vectors for hightemperature T > Tc simulations and 128 noise vector for low temperature T < Tc
simulations. An example calculation is shown as follows. See Fig. 27. This figure
shows that the noise method can produce consistent results with these obtained by LU
decomposition within a range of statistical errors. We thus use the noise method for the
trace of winding number expansion.
Note that we can reduce the computation cost of WNE using this method. This isbecause the times of loops are reduced as Nx Ny Nz Nt × 3 × 4 = 24576 → Nnoise = 128
whensize is) 83 × 4. Then, the time complexity of WNE with noise trace is( the lattice2
O (Nx Ny Nz ) × (Nt )2 .

D∗

Direct method via winding number expansion without
canonical ensembles

When we calculate thermodynamic observables using the canonical approach with WNE,
we calculate the following quantities in following order,
Wk → Z(iµI ) → ZC (n) → Z(µ) .
53

(133)

Figure 27: Pure imaginary chemical potential dependence of fermion determinants. Red
and blue points are calculated by the winding number expansion with 16 noise vectors
and LU decomposition respectively.
Then we use WNE, Fourier transform and fugacity expansion. We can avoid the sign
problem using above method.
On the other hand, we notice that we can calculate Z(µ) from Wk directly,
Wk → Z(µ) ,(134)

because we can use the WNE of fermion determinant at any µ (not only iµI ). (SeeFig. 28.) In this thesis, we call it the direct method.
It may seem that the direct method works out well without sign problem. However,this method diﬀers from the standard canonical approach in the following points.
• In the direct method, we calculate Z(µ) without the canonical ensembles unlike
the canonical approach.
• In the direct method, we perform the finite density calculation,[]∑
det ∆(µ) = exp
Wk (eµ/T )k ,

(135)

k

before we calculate the integral of gauge field. On the other∑hand, in the canonicalapproach, we perform the finite density calculation Z(µ) = n ZC (n)(eµ/T )n after
we calculate the integral.
From the second point, we expect that the direct method
does not work well because the∫sign problem occurs when we calculate Z(µ) = DU [det ∆(µ)]Nf eSG using complex
det ∆(µ).
54Figure 28: Diﬀerence of canonical approach and direct method.
Actually, the canonical approach is more accurate than the direct method. For
example, although the pressure is real, sometimes it has the imaginary part since there
are numerical errors. Figure 29 shows the imaginary part of pressure which we calculate
using the canonical approach and direct method. From this figure, we see that the resultsof canonical approach and direct method are consistent to zero. However, the statisticalerror of canonical approach is smaller than the one of direct method. We can see thisproperty from other observables. For instance, Fig. 30 and 31 show the imaginary partof baryon number density and baryon susceptibility. We find that the error of canonicalapproach is always smaller than the one of direct method.
Figure 32, 33 and 34 show the thermodynamic observables which we calculate using the canonical approach and direct method. Then, the validity range of canonicalapproach is estimated as Sec. 8.3.
The range of direct method is estimated by the magnitude of statistical error ofimaginary part. The procedure is as follows.

55

Figure 29: Imaginary part of pressure of canonical approach and direct method. The red
and green points are the data of direct method and canonical approach at T /Tc = 1.08,
respectively.

Figure 30: Imaginary part of baryon number density of canonical approach and direct
method. The red and green points are the data of direct method and canonical approachat T /Tc = 1.08, respectively.

56

Figure 31: Imaginary part of baryon susceptibility of canonical approach and direct
method. The red and green points are the data of direct method and canonical approachat T /Tc = 1.08, respectively.

10
Canonical Direct

8
6
4
2
0

0

1

2

3

4

5

Figure 32: Baryon chemical potential dependence of pressure which is calculated bycanonical approach and direct method at T /Tc = 1.35, 1.20, 1.08, 0.99, 0.93 and 0.84.

57

3.0

Canonical Direct

2.5
2.0
1.5
1.0
0.5
0.0

0.0

1.0

2.0

3.0

4.0

Figure 33: Baryon chemical potential dependence of baryon number density which is
calculated by canonical approach and direct method at T /Tc = 1.35, 1.20, 1.08, 0.99,
0.93 and 0.84.

1.2

Canonical Direct

1.0
0.8
0.6
0.4
0.2
0.0

0.0

1.0

2.0

3.0

Figure 34: Baryon chemical potential dependence of baryon susceptibility which is calculated by canonical approach and direct method at T /Tc = 1.35, 1.20, 1.08, 0.99, 0.93
and 0.84.

58

1. Save the error of imaginary part
observables of canonicalapproach at the bound( of[] )4ary of validity range, e.g. δ Im ∆p(µval(can) )/T can .
( [] )2. Compare the one of direct method δ Im ∆p(µ)/T 4 dir with)( []
δ Im ∆p(µval(can) )/T 4 can .
3. Find µ1 which satisfies,
( [
( [
] )
] )
δ Im ∆p(µ1 )/T 4 dir = δ Im ∆p(µval(can) )/T 4 can .

(136)

4. Consider the boundary of validity range of direct method is µval(dir) = µ1 .
That is, we plot the data of canonical approach and direct method with the same accuracy.
From these figure, we find that the validity range of canonical approach is larger thanthe one of direct method. Especially, in the pressure data, canonical approach is valid
beyond µB /T = 5 at T /Tc = 1.20, although direct method is valid to µB /T = 4.
Therefore, when we calculate the thermodynamic observables at finite density, it isbetter that we use the canonical approach in order to avoid the sign problem.

E∗

Algorithm for estimation of cancellation of significant
digitsIn this section, we will see how to estimate the cancellation of significant digits. The
cancellation is caused by the subtraction between two similar numbers.
The logarithm of numbers give us the magnitude of the numbers. For example, thelogarithm of a = 6.022 × 1023 is log10 a = 23.78. Then, “23” corresponds to the power ofa, and “0.78” corresponds to the significant digits of a. This is because we can dissolve
log10 a as
log10 a = log10 (6.022 × 1023 ) = log10 (6.022) + log10 (1023 ) .
(137)
Using this idea, we can estimate how many digits are the same between two numbers.
Let a and b be real numbers. For simplicity, we consider that a and b are positive
and a > b. Here, we define the following quantity,
r := a − b ,(138)

d := log10 (a) − log10 (r) .

(139)

where “r” means the result and “d” means the diﬀerence between the magnitudes of aand r. If a and b is
a = 6.022 × 1023 ,(140)

b = 6.021 × 10

(141)

59

23

,

then r and d becomes the following number,r = log10 (0.001 × 1023 ) = 20.00 ,(142)

d = 23.78 − 20.00 = 3.78 .

(143)

In this case, when we calculate r = a − b, three significant digits are lost. This
“three” appears in Eq. (143). That is, the number of canceled digits are given as(Number of canceled digits) = ⌊d⌋ = ⌊log10 (a) − log10 (r)⌋ ,(144)

where ⌊·⌋ is the floor function.
Here, what means “0.78” of Eq. (143)? This is the diﬀerence between the magnitudes of significant digits. Let asig and apow be the significant digits and power of a,respectively. We then can write d as
d = log10 (a) − log10 (r) = log10 (asig × 10apow ) − log10 (rsig × 10rpow ))(( apow )asig
10+ log10= log10rsig
10rpow()
asig= log10+ (apow − rpow ) ,rsig⌊d⌋ = apow − rpow .

(145)
(146)

Thus, “0.78” is log10 (asig /rsig ). (In fact, we can compute it as log10 (6.022/ 1.0) = 0.78.)
Note that 1 ≤ asig /rsig < 10. This is because asig and rsig satisfies 1 ≤ asig < b and1 ≤ rsig < b where b is the base of logarithm. (b = 10; See also Sec. 7.3.)
The Algorithm for estimation of the number of canceled digits is as follows. Notethat we then adopt the multi precision calculation.
1. Preparation: adopt the multi precision calculation to the program.
2. Initialization:(a) Set the number of significant digits is large enough.
(b) Choose a as the minuend and b as subtrahend. They are multi precisionvariables.
3. Estimation:(a) Save log10 (a).
(b) Calculate the subtraction r = a − b which may cause the cancellation ofsignificant digits.
(c) Calculate d = log10 (a) − log10 (r).
(d) Print ⌊d⌋ which is the number of canceled digits.

60

References[1] K. Fukushima and T. Hatsuda, “The phase diagram of dense QCD,” Rept. Prog.
Phys. 74 (2011) 014001.
[2] M. Buballa, and M. Oertel, “Color–flavor unlocking and phase diagram with self–consistently determined strange-quark masses,” Nucl. Phys. A 703 (2002) 770.
[3] C. Itzykson, and J. B. Zuber, “Quantum Field Theory,” Dover Publications (2006).
[4] E.D. Bloom, D.H. Coward, H.C. DeStaebler, J. Drees, G. Miller, L.W. Mo, R.E.
Taylor, M. Breidenbach, J.I. Friedman, G.C. Hartmann and H.W. Kendall, “HighEnergy Inelastic e–p Scattering at 6◦ and 10◦ ,” Phys. Rev. Lett. 23 (1969) 930.
[5] J.D. Bjorken and E.A. Paschos, “Inelastic Electron–Proton and γ–Proton Scatteringand the Structure of the Nucleon,” Phys. Rev. 185 (1969) 1975.
[6] D.J. Gross, F. Wilczek. “Ultraviolet behavior of non-abelian gauge theories,” Phys.
Rev. Lett. 30 (1973) 1343.
[7] H.D. Politzer. “Reliable perturbative results for strong interactions,” Phys. Rev.
Lett. 30 (1973) 1346.
[8] M. Creutz, “Monte Carlo study of quantized SU(2) gauge theory,” Phys. Rev. D 21
(1980) 2308.
[9] M. Creutz, “Asymptotic-freedom scales,” Phys. Rev. Lett. 45 (1980) 313.
[10] D.J. Gross, R.D. Pisarski and L. G. Yaﬀe, “QCD and instantons at finite temperature,” Rev. Mod. Phys. 53 (1981) 43.
[11] A. Roberge and N. Weiss. “Gauge theories with imaginary chemical potential andthe phases of QCD,” Nucl. Phys. B 275 (1986) 734.
[12] Z. Fodor and S.D. Katz, “Critical point of QCD at finite T and μ, lattice resultsfor physical quark masses,” JHEP 0404 (2004) 050.
[13] S. Datta, R.V. Gava and, S. Gupta, “The QCD Critical Point : marching towards
continuum,” Nucl. Phys. A 904 (2013) 883.
[14] O. Scavenius, A. Mocsy, I.N. Mishustin and D.H. Rischke, “Chiral Phase Transition
within Eﬀective Models with Constituent Quarks,” Phys. Rev. C 64 (2001) 045202.
[15] K. Jansen and R. Sommer, “O(a) improvement of lattice QCD with two flavors of
Wilson quarks,” Nucl. Phys. B 530 (1998) 185.
[16] Y. Iwasaki, “Renormalization Group Analysis of Lattice Theories and Improved Lattice Action. II — four–dimensional non–abelian SU(N) gauge model —,” UTHEP118 (1983), preprint.
61

[17] B. Sheikholeslami and R. Wohlert, “Improved continuum limit lattice action for
QCD with Wilson fermions,” Nucl. Phys. B259 (1985) 572.
[18] A. Nakamura, “Quarks and Gluons at Finite Temperature and Density,” Phys. Lett.
149B (1984) 391.
[19] Ph. de Forcrand, “Simulating QCD at finite density,” PoS (LATTICE2009) 010(2010).
[20] J. Engels and H. Satz, “Deconfinement at finite baryon number density,” Phys.
Lett. B 159 (1985) 151.
[21] I. Barbour, N.E. Behilil, E. Dagotto, F. Karsch, A. Moreo, M. Stone, and H.W.
Wyld, “Problems with finite density simulations of lattice QCD,” Nucl. Phys. B
275 (1986) 296.
[22] Z. Fodor and S.D. Katz, “A new method to study lattice QCD at finite temperatureand chemical potential,” Phys. Lett. B 534 (2002) 87.
[23] F. Csikor, G.I. Egri, Z. Fodor, S.D. Katz, K.K. Szabo and A.I. Toth, “The QCD
equation of state at finite T/μ on the lattice,” Prog. Theor. (Phys. Suppl.) 153(2004) 93.
[24] K. Nagata and A. Nakamura, “EoS of finite density QCD with Wilson fermions by
Multi–Parameter Reweighting and Taylor expansion,” JHEP 04 (2012) 092.
[25] S. Ejiri, C.R. Allton, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and C.
Schmidt, “Study of QCD thermodynamics at finite density by Taylor expansion,”
Prog. Theor. (Phys. Suppl.) 153 (2004) 118.
[26] A. Bazavov, H.T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y.
Maezawa, S. Mukherjee, H. Ohno, P. Petreczky, H. Sandmeyer, P. Steinbrecher, C.
Schmidt, S. Sharma, W. Soeldner and M. Wagner, “The QCD Equation of State toO(µ6B ) from Lattice QCD,” Phys. Rev. D 95 (2017) 054504.
[27] C. Schmidt and S. Sharma, “The phase structure of QCD,” J. Phys. G 44 (2017)
104002.
[28] WHOT–QCD Collaboration: S. Ejiri, Y. Maezawa, N. Ukita, S. Aoki, T. Hatsuda,
N. Ishii, K. Kanaya and T. Umeda, “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 (2010) 014508.
[29] A. Hasenfratz and D. Toussaint, “Canonical ensembles and nonzero density quantum chromodynamics,” Nucl. Phys. B 371 (1992) 539.
[30] M.M. Aggarwal, et al, “Higher Moments of Net Proton Multiplicity Distributions
at RHIC,” Phys. Rev. Lett. 105 (2010) 022302.
62

[31] A. Nakamura, R. Fukuda, S. Oka, S. Sakai, A. Suzuki, Y. Taniguchi, K. Morita and
K. Nagata, “Net-Baryon Multiplicity and QCD Phase Diagram,” PoS (CPOD2014)
021 (2015).
[32] Y. Sasai, A. Nakamura and T. Takaishi, “Phase Fluctuation of Fermion Determinant
in Lattice QCD at Finite Density,” Nucl. Phys. B (Proc. Suppl.) 129–130 (2004)
539[33] S. Kratochvila and Ph. de Forcrand, “The canonical approach to Finite Density
QCD,” PoS (LATTICE2005) 167 (2006).
[34] Ph. de Forcrand and S. Kratochvila, “Finite density QCD with a canonical approach,” Nucl. Phys. B (Proc. Suppl.) 153 (2006) 62.
[35] X. Meng, A. Li, A. Alexandru and K.F. Liu, “Winding number expansion for thecanonical approach to finite density simulations,” PoS (LATTICE2008) 032 (2009).
[36] 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 (2012) 014502.
[37] K. Nagata and A. Nakamura, “Wilson Fermion Determinant in Lattice QCD,” Phys.
Rev. D 82 (2010) 094027.
[38] A. Alexandru and U. Wenger, “QCD at non–zero density and canonical partition
functions with Wilson fermions,” Phys. Rev. D 83 (2011) 034502.
[39] P.E. Gibbs, “The fermion propagator matrix in lattice QCD,” Phys. Lett. B 172(1986) 53.
[40] J. Danzer and C. Gattringer, “Winding expansion techniques for lattice QCD withchemical potential,” Phys. Rev. D 78 (2008) 114506.
[41] A. Nakamura, S. Oka and Y. Taniguchi, “QCD phase transition at real chemical
potential with canonical approach,” JHEP 02 (2016) 054.
[42] R. Fukuda, A. Nakamura, S. Oka, S. Sakai, A. Suzuki and Y. Taniguchi, “Beating
the sign problem in finite density lattice QCD,” PoS (LATTICE2015) 208 (2016).
[43] H. Nyquist, “Certain topics in telegraph transmission theory,” Trans. AIEE. 47(1928) 617.
[44] C. E. Shannon, “Communication in the Presence of Noise,” Proceedings of the IRE
37 (1949) 10.
[45] S. Oka, “Exploring finite density QCD phase transition with canonical approach—Power of multiple precision computation—,” PoS (LATTICE2015) 067 (2016).

63

[46] R. Fukuda, A. Nakamura and S. Oka, “Canonical approach to finite density QCD
with multiple precision computation,” Phys. Rev. D 93, 094508 (2016).
[47] D. E. Knuth, “The Art of Computer Programming, Volume 2 / Seminumerical
Algorithms,” Second Edition, Addison–Wesley publishing (1969).
[48] D. Goldberg, “What every computer scientist should know about floating-point
arithmetic,” ACM Computing Surveys (CSUR) 23 (1991) 5.
[49] A.A. Khan et al. (CP-PACS Collaboration), “Equation of state in finite temperature
QCD with two flavors of improved Wilson quarks,” Phys. Rev. D 64 (2001) 074510.
[50] R. Fukuda, A. Nakamura and S. Oka, “Validity range of canonical approach to finite
density QCD,” PoS (LATTICE2015) 167 (2016).
[51] C. Gattringer and T. Kloiber, “Lattice study of the Silver Blaze phenomenon for acharged scalar ϕ4 field,” Nucl. Phys. B 869 (2013) 56.
[52] Ph. de Forcrand and O. Philipsen, “The QCD phase diagram for small densities
from imaginary chemical potential,” Nucl. Phys. B 642 (2011) 290.
[53] R. Fukuda, A. Nakamura, S. Oka, A. Suzuki and Y. Taniguchi, “Study of the sign
problem in canonical approach,” PoS (LATTICE2016) 060 (2017).
[54] V.A. Goy, V. Bornyakov, D. Boyda, A. Molochkov, A. Nakamura, A. Nikolaev andV. Zakharov, “Sign problem in finite density lattice QCD,” PTEP 2017 (2017)
031D01.
[55] A. Nakamura, S. Oka, A. Suzuki and Y. Taniguchi, “Calculation of high–order
cumulants with canonical ensemble method in lattice QCD,” PoS (LATTICE2015)
168 (2016).
[56] D.L. Boyda, V.G. Bornyakov, V.A. Goy, V.I. Zakharov, A.V. Molochkov, A. Nakamura and A.A. Nikolaev, “Novel approach to deriving the canonical generating
functional in lattice QCD at a finite chemical potential,” JETP Lett. 104 (2016)657.
[57] 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 Nf = 2 lattice QCD at finite baryon density,” Phys. Rev. D 95 (2017)094506.
[58] D. Boyda, V.G. Bornyakov, V. Goy, A. Molochkov, A. Nakamura, A. Nikolaev and
V.I. Zakharov, “Lattice Study of QCD Phase Structure by Canonical Approach,”
arXiv:1704.03980 [hep-lat], preprint.
[59] A. Alexandru, M. Faber, I. Horvath and K.F. Liu, “Lattice QCD at finite density
via a new canonical approach,” Phys. Rev. D 72 (2005) 114513.
64[60] A. Li, A. Alexandru, K.F. Liu and X. Meng, “Finite density phase transition of
QCD with Nf = 4 and Nf = 2 using canonical ensemble method,” Phys. Rev. D
82 (2010) 054502.
[61] A. Li, A. Alexandru and K.F. Liu, “Critical point of Nf = 3 QCD from latticesimulations in the canonical ensemble,” Phys. Rev. D 84 (2011) 071503.
[62] A. Alexandru, C. Gattringer, H.P. Schadler, K. Splittorﬀ and J.J. M. Verbaarschot,“Distribution of canonical determinants in QCD,” Phys. Rev. D 91 (2015) 074501.
[63] C. Gattringer, H.P. Schadler, “Generalized quark number susceptibilities from fugacity expansion at finite chemical potential for Nf = 2 Wilson fermions,” Phys.
Rev. D 91 (2015) 074511.
[64] T. Matthews and A. Salam, “Propagators of Quantized Field,” Nuovo Cimento 2(1955) 120.
[65] J. Foley, K.J. Juge, A. O’Cais, M. Peardon, S.M. Ryan and J.I. Skullerud, “Practical
all–to–all propagators for lattice QCD,” Comput. Phys. Commun. 172 (2005) 145.

65