Canonical approach -- Investigation of finite density QCD phase transition
Shotaro Oka
hep-lat physics.comp-ph





















































empty

left=0.5cm, right=0.5cm
2.6cm



	Canonical approach

	


	— Investigation of finite density QCD phase transition —





	


	Shotaro Oka
	
	

	
	Rikkyo University




	


	A thesis submitted for the degree of

	

	
	Doctor of Philosophy

	


	January 2018













empty

-2.6cm
0.3cm
0.3cm

































QCD is the fundamental theory which describes the dynamics of quarks and gluons.
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 the studies 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 approach
of 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 to enhance the calculation accuracy of physical quantities,
we adopt the multi precision calculation to our program based on the canonical approach.
In this thesis, we will 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 do not see the QCD phase transition yet.
However, in this thesis, we find that canonical approach 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 window of study of QCD phase diagram at finite density.
With our improvement, canonical approach has the possibility for investigation of thermodynamic observables
at any chemical potential.

  













* 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 to finite density QCD,"
		PoS (LATTICE2015) 167 (2016).








plain





§ INTRODUCTION






 §.§ 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.




However, we have expected that QCD has a rich phase structure
at finite temperature and 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 made of hadronic matters.

At high temperature, T > 10^2 MeV (= 10^12 K)<cit.>, quarks are disentangled from hadrons.
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 the 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 the early universe for ages.

At high baryon density, μ_B > 10^3 MeV (= 10^18 kg/m^3)<cit.>,


quarks can not distinguish the boundary of baryons (see Fig. <ref>).
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.





We have mainly obtained a lot of interesting information about the quark matter
by 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.









 §.§ 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.
As an electric charge is origin of the electromagnetic interaction,
color charge causes the strong 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

    A_μ = A^a b_μ(x,y,z,t) = A^α_μ(x,y,z,t) T^α, a b  .


Then, QCD Lagrangian is written as follows,

    L    = -1/4 G_μν^α G^αμν + ψ̅^a (i γ^μ D_μ - m) ψ^a    ,  
    
    	G^α_μν   = ∂_μ A^α_ν - ∂_ν A^α_μ + g f^αβγ A^β_μ A^γ_ν  ,  
    
    	D_μψ^a    = ∂_μψ^a - i g A^α_μ T^α, a bψ^b    ,

where G^α_μν is the strength of the gauge field,
D_μ is the covariant derivative and g is the coupling 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^β]    = i f^αβγ T^γ  .


Here, since gluons have color indices (that is, gluons are not commutative), it can interact with itself.
We can see this from the first term of Eq. (<ref>),

    G_μν^α G^αμν   = ( ∂_μ A^α_ν - ∂_ν A^α_μ + g f^αβγ A^β_μ A^γ_ν)  
       ×( ∂^μ A^α ν - ∂^ν A^α μ + g f^αδϵ A^δ μ A^ϵ ν)  
       = ∂_μ A^α_ν∂^μ A^α ν
    			- ∂_μ A^α_ν∂^ν A^α μ
    			+ g ∂_μ A^α_ν f^αδϵ A^δ μ A^ϵ ν
       = - ∂_ν A^α_μ∂^μ A^α ν
    			+ ∂_ν A^α_μ∂^ν A^α μ
    			- g ∂_ν A^α_μ f^αδϵ A^δ μ A^ϵ ν
       = + g f^αβγ A^β_μ A^γ_ν∂^μ A^α ν
    			- g f^αβγ A^β_μ A^γ_ν∂^ν A^α μ
       = + g^2 f^αβγ A^β_μ A^γ_ν f^αδϵ A^δ μ A^ϵ ν  .

The underlined terms of Eq. (<ref>) describe the interactions between gluons,
and the other terms 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 affects the dynamics of quarks significantly.
For example, quarks (color charges) emit or absorb the gluons
as well as electric charges do for the photons.
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 such a 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 in the shortest way,
and a quark can not appear alone since it needs infinite energy to be free.

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. <ref>; it is a meson case.)

This phenomenon is called the confinement of quarks.




QCD also has other interesting properties,
for instance, the asymptotic freedom is a significant 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,

    B(g(μ))    = μd g(μ)/d μ = -B_0 g^3 - B_1 g^5 + ⋯  ,  
    
    	B_0    = 1/(4 π)^2( 11/3 N_c - 2/3 N_f )    ,  
    
    	B_1    = 1/(4 π)^2( 34/3N_c^2 - 10/3 N_c N_f - N_c^2 - 1/N_c N_f )    ,

using the two–loop perturbation theory<cit.>.
Here, N_c = 3 is the number of colors and N_f is the number of flavors.
From these, we find the following relation,

    log(μ) = 1/2 B_0 g^2 + B_1/2 B_0^2log(g^2) - B_1/2 B_0log(B_0 + B_1 g^2) + (Const.)    .

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,
and asks us what happens between hadron gas phase and QGP phase;
these two phases are respectively ruled by different particles and physics.

The asymptotic freedom was found by a collider experiment of SLAC<cit.> in 1968,
which was explained theoretically by D.J. Gross and F. Wilczek<cit.>,
and H.D. Politzer<cit.> in 1973.



In addition, it is confirmed that the results of string tension[
	We can consider the gluon flux tube which bonds two quarks as a string.
	Then, introducing the string 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.
] of a lattice QCD simulation
reproduce the theoretical expectation with high accuracy as Fig. <ref>.
That is, the non–perturbative result supports the perturbative exception
and shows the behavior of running coupling constant at middle region
where the perturbation theory does not work;
it is a valuable result.
This simulation was performed by M. Creutz<cit.> in 1980.















 §.§ Finite temperature and density QCD


When we discuss actual QCD systems, we use methods of statistical mechanics and thermodynamics.
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 the following way.

In order to introduce temperature to QCD,
we consider a partition function,

    Z = Tr( e^-β ( Ĥ - μN̂ ))    .

where β = 1/T is an inverse temperature, Ĥ is the Hamiltonian operator
and N̂ is the number operator.
Let | ϕ⟩ be a state of the system on the Hilbert space.
Then, we can rewrite the trace as a summation of physical states,

    Z = ∑_ϕ⟨ϕ(x) | e^-β ( Ĥ - μN̂ )| ϕ(x) ⟩  .

For simplicity, we will neglect fermions and μN̂ for a moment.
Then, we notice that the operator of Eq. (<ref>), exp(-βĤ),
corresponds to the Euclidean time evolution operator of quantum mechanics, exp(-Ĥ t).
We therefore can consider that, essentially, β is the elapsed time of states.

Now, states are represented by gauge field A_μ.
Taking the Coulomb gauge, we can write the partition function as

    Z = ∑_𝐀⟨𝐀(β, 𝐱) | e^-βĤ| 𝐀(0, 𝐱) ⟩  ,

where 3–dimensional vector 𝐀 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. (<ref>).
In order to remove this inconsistency, we require the periodic boundary condition along time direction,

    𝐀(β, 𝐱) = 𝐀(0, 𝐱)    .

Requiring this condition, we can construct the finite temperature QCD.
(See also Ref. <cit.>.)



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. (<ref>).
Using the path integral formalism, we can write the partition function
with chemical potential μ as follows,

    Z = ∫ DA Dψ Dψ̅exp( -∫_0^β dτ∫ d^3x [ -1/4 F^2 + ψ̅ (γ D - m) ψ - μψ^†ψ] )   ,

where τ = i t 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̂ = ∫ d^3x ψ^†ψ
and N̂ is constant of time.






































































However, it is complicated because ψ^† appears.
We can simplify this equation via the change of variables

    A_ν→ A_ν - i/g( ∂_ν e^μτ) e^-μτ  .

The partition function now becomes

    Z = ∫ DA Dψ Dψ̅exp( -∫_0^β dτ∫ d^3x [ -1/4 F^2 + ψ̅ (γ D - m) ψ] )    .







That is, the chemical potential appears as the imaginary part of temporal component of gauge field.

A. Roberge and N. Weiss considered above discussion
with pure imaginary chemical potential μ = i μ_I<cit.>.
Then, because e^i μ_I is periodic, Z(μ_I) has the periodicity,

    Z(μ_I) = Z( μ_I + 2 π k/N)    ,

for any integer k.
This result is proved for SU(N) theory.
(N=3 is QCD case.)
Thus, the partition function Z(μ) of QCD has Z_3 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 numerical methods.
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 T_c ≈ 150 MeV
from zero density lattice QCD simulations[
	The critical point is located at finite density.
	This means that we can not use usual order parameters such as the Polyakov loop and chiral condensate
	because the center symmetry of SU(3) and chiral symmetry 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 T_pc which is approximately estimated at finite density
	by the Polyakov loop or chiral condensate.
	In this thesis, we write the both of critical and pseudo critical temperature as T_c.
]
<cit.>.
But, critical density (chemical potential) μ^c_B is not known well.

Z. Fodor and S.D. Katz studied it in 2004 for 2+1 flavor lattice QCD<cit.>,
then they gave μ^c_B = 320(40) MeV.
On the other hand, S. Datta, R.V. Gavai and S. Gupta studied it in 2012 for two flavor lattice QCD<cit.>.
They then gave μ^c_B = 243(20) MeV.
Moreover, O. Scavenius, A. Mocsy, I.N. Mishustin and D.H. Rischke introduced
the phase diagrams of σ–model and NJL model which are effective models of QCD, in 2000<cit.>.
From this paper, critical point of these models are at μ^c_B = 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 (and several models) can not estimate physical quantities well
in this region of finite chemical potential.

Although the sign problem is a severe difficulty,
ones have been studying how to avoid this problem.
As one of such methods, the canonical approach was developed.
This method 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 effect of the sign problem
and investigate the QCD phase diagram at finite density by solving the problems of canonical approach.









§ BASIC DEFINITIONS OF FINITE DENSITY LATTICE QCD





 §.§ 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, fermions are on the sites of the lattice,
and gluons are on the links of the lattice because they have a direction as Lorentz index.


Let a be the lattice spacing.
Because A_μ is not covariant under the SU(3) gauge transformation,
we alternatively use the following covariant quantity,

    U_μ (x) = U (x → x+μ̂a) := e^i a g A_μ(x+μ̂/2)∈SU(3)  ,

for a lattice variable.
Here, μ̂ is the unit vector along μ direction.
U_μ are called the link variables.

Since they are unitary matrices, we find the following relation,

    U_μ(x) U_-μ(x+μ̂a) = 1 = U_μ(x) U^†_μ(x)    ,

thus U_-μ(x+μ̂a) = U^†_μ(x).




We can write the gauge action by gauge invariant quantities which are constructed by U_μ.
For example, the trace of closed loops of U_μ is gauge invariant.
The simplest loop, that is a unit square on the lattice,

    tr U_μ(x) U_ν(x+μ̂) U_μ^†(x+ν̂) U_ν^†(x)    ,

is called the plaquette.

It is invariant under the following gauge transformation,

    Ψ(x)    →Ψ^'(x) = e^-i θ^α(x) T^αΨ(x) =: Ω(x) Ψ(x)    ,  
    
    	A_μ(x)    → A^'_μ(x) = Ω(x) A^α_μ T^αΩ^†(x) - i/g [∂_μΩ(x)] Ω^†(x)    ,

with Ω(x) ∈ SU(3). That is,

    U_μ(x) → U^'_μ(x) = Ω(x) U_μ(x) Ω^†(x+μ̂a)    .

Using the plaquette, we can write the gauge action as

    S_G = β/3∑_μ≠νtr U_μ(x) U_ν(x+μ̂) U_μ^†(x+ν̂) U_ν^†(x)    .

Now, we use β as the effective coupling constant. Note that it is not the inverse temperature.
This action is called Wilson gauge action.


It is difficult to construct the fermion action.
This is because fermions get extra degrees of freedom
when we discretize the Dirac operator naively.
This problem is known as the fermion doubling problem.

To reduce the problem, K.G. Wilson considered the following action,

    S_F    = a^4 ∑_x, x^'Ψ̅(x) Δ(x, x^') Ψ(x^')    ,  
    Δ(x, x^')    = δ_x, x^'
    		- κ∑_i=1^4( ( r - γ_i ) U_i (x) δ_x^', x+î
                                         + ( r + γ_i ) U_i^† (x^') δ_x^', x-î)     ,

where κ is the hopping parameter.
κ rules the movement of quarks, which depends on the quark mass m_q<cit.> as

    κ = 1/8 + 2 m_q 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.







 §.§ 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^-S_G e^-S_F(μ) .

Then U is link variables, Ψ and Ψ̅ are fermions and anti–fermions
on a lattice respectively.
S_G and S_F are a gauge action and a fermion action respectively.
We will see them in detail below. 




For the gauge part S_G, we use the improved Wilson gauge action,

    S_G = β/3∑_μ≠ν[ c_0 tr P_μν^1 × 1 + c_1 tr P_μν^1 × 2]   ,

where β = 6/g^2 is the effective coupling constant (g is a bare coupling constant),
P_μν^1 × 1 is the usual plaquette
and P_μν^1 × 2 is the rectangle plaquette:

    P_μν^1 × 1 (x)    = U_μ(x) U_ν(x+μ̂) U_μ^†(x+ν̂) U_ν^†(x)  
    
    	P_μν^1 × 2 (x)    = U_μ(x) U_μ(x+μ̂) U_ν(x+2μ̂) 
         U_μ^†(x+μ̂+ν̂) U_μ^†(x+ν̂) U_ν^†(x)  .

We set the coefficients in Eq. (<ref>) as c_1 = -0.331 and c_0 = 1 - 8 c_1 as Ref. <cit.>.
It is called Iwasaki gauge action. 


This action was suggested by Y. Iwasaki in 1983.
The aim of this action is that
the numerical 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 S_F is written as 

    S_F(μ) = a^4 ∑_x, x^'Ψ̅(x) Δ(μ, x, x^') Ψ(x^')

where Δ(μ) is a ferimon matrix.
Now, we adopt the improved Wilson fermions with chemical potential μ for future reference,

    Δ(μ, x, x^')    = δ_x, x^'
    		- κ∑_i=1^3( ( 1 - γ_i ) U_i (x) δ_x^', x+î
                                                 + ( 1 + γ_i ) U_i^† (x^') δ_x^', x-î)  
         - κ( e^+μ a( 1 - γ_4 ) U_4 (x) δ_x^', x+4̂
                                  + e^-μ a( 1 + γ_4 ) U_4^† (x^') δ_x^', x-4̂)  
         - κ C_SWδ_x, x^'∑_μ < νσ_μν F_μν (x)   ,
  
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),

    F_μν = P_μν^1 × 1 (x) + P_ν -μ^1 × 1 (x) + P_-μ -ν^1 × 1 (x) + P_-νμ^1 × 1 (x)    ,

and σ_μν = i [γ_μ, γ_ν] / 2.


Then, we introduce e^μ a to the terms of U_4 and U^†_4
because of Eq. (<ref>). 




This means that fermions get the coefficient e^μ a when they travel
along t–coordinate.
Since the inverse temperature is related with the boundary conditions of ψ and A_μ,
if fermions wind around the lattice along the time direction,
they get e^μ a N_t = e^μ / T.


The term of C_SW 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 the discretization error of the action.
By a result of one–loop perturbation theory<cit.>,
we set C_SW = (1 - 0.8412/β)^-3/4. (This β is the coupling constant.)







Note that we can integrate out e^-S_F

because it forms the Gaussian integral of the Grassmann number Ψ.
(See Appx. <ref>.)
As a result, Eq. (<ref>) is rewritten as follows.

    Z(μ) = ∫ DU [ Δ(μ) ] e^-S_G  .

It is obtained for one flavor fermion case.
If we want to consider many flavor fermions,
we can use the following equation, 

    Z(μ) = ∫ DU [ Δ(μ) ]^N_f e^-S_G  ,

where N_f is the number of fermions. 








Using above definitions, we can calculate the expectation values of thermodynamic observables,

    ⟨ O(μ) ⟩ = 1/Z(μ)∫ DU O(μ) [ Δ(μ) ]^N_f e^-S_G .

Here, operater O(μ), fermion matrix Δ(μ) and gauge action S_G
depend on link variables U.
Therefore, in a practical caluculation,
we calculate them for each configuration (a set of { U }) and average their product.









§ BASIC METHODS OF FINITE DENSITY LATTICE QCD





 §.§ Monte Carlo method



Here, link variables U have many indices, specifically

    U = U^a b_μ(x,y,z,t)    ,

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 × 10^5 times integration.



Moreover, we have to perform ∫ D U^a b_μ(x,y,z,t) using the division quadrature.
That is,
if we write the component of gauge field with fixed a, b, μ and (x,y,z,t)
as U^a b_μ(x,y,z,t) =: u,

    ∫ d (U^a b_μ(x,y,z,t)) f[U]    = ∫ d u f(u)  
       ≈∑_k=1^N f(u_k) Δ u    ,

where u_k (k ∈{ 1, 2, ⋯ N }) is a sampled u.
Therefore, if N=10, we should perform 10^3.6 × 10^5 operations in this case.
Hence, it is difficult to calculate Eq. (<ref>) naively. 




We then can use the Monte Carlo method in order to reduce
the computational cost of 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 as 

    p ∝[ Δ(μ) ]^N_f e^-S_G  ,

in this case.
We are going to see it in detail below.





 §.§ Improved Monte Carlo methods






  §.§.§ 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 efficient method because it generates
not only statistically important configurations
but also unimportant others 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 it as a filter.
For instance, we consider a simple Monte Carlo case,

    E[O] = ∫ O(x) f(x) dx    ,

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,

    E[O] = lim_m →∞1/m∑_j=1^m O(x_m)

where x_m 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. (<ref>) as

    E[O] = ∫ O(x) f(x)/f^'(x) f^'(x) dx    .

Thus, Eq. (<ref>) can be rewritten,

    E[O] = lim_m →∞1/m∑_j=1^m O(x^'_m) f(x^'_m)/f^'(x^'_m)  ,

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.





  §.§.§ 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 the Metropolis 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)    .


Metropolis algorithm is constructed as follows.

	
  * Initialization: choose an arbitrary point x_0 as an initial sample and
		an arbitrary probability density Q(x|y).

	
  * Iteration: for each t,
	
		
  * Generate new sample x^' from Q(x^'|x_t).

		
  * Calculate acceptance ratio α = f(x^')/f(x_t).
			Note that it is just the weight function.
		
		
  * If α≥ 1, x^' is more likely than x_t,
			then it is accepted for next sample x_t+1;
			if α < 1, x^' is accepted according to probability α;
			when x^' is rejected, set x_t+1 = x_t.
	


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, this algorithm is not efficient.
Usually, an improved method is alternatively used;
we will see it in the following section.




  §.§.§ 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.




	
  * Initialization:
	
		
  * Choose ϕ as the integration variable.
		
  * Generate the conjugation momentum π with
			probability P(π) ∝ e^-π^2/2 (Gauss distribution).
	
	
	
  * Iterations (HMD part):
		evolve ϕ and π through the following equations,
		
    {[        H[ϕ, π] = -1/2π^2 + S(ϕ);           ϕ̇ = ∂ H[ϕ, π]/∂π = π; π̇ = -∂ H[ϕ, π]/∂ϕ = -∂ S(ϕ)/∂ϕ ].

		where H[ϕ, π] is the Hamiltonian and
		S(ϕ) is the action for ϕ (it is given).

	
  * Metropolis test: using the Metropolis method,
		accept or reject ϕ and π according to the following probability,
		
    P( {ϕ_init, π_init}→{ϕ_new, π_new} ) = min{ 1, e^-Δ H}  ,  
       Δ H = H(ϕ_new, π_new) - H(ϕ_init, π_init)    .


	
  * Repeat 2 and 3.






















§ SIGN PROBLEM







 §.§ Sign problem


Basically, using above definitions and methods,
we can start the calculation of finite density lattice QCD.
However, Δ becomes a complex number when μ is nonzero.
We can see this using the definition of Δ and the property of the determinant.
First, we observe the following relation by a property of a determinant:

    [ Δ(μ) ]^⋆ = Δ^† (μ)   .

Then, Δ^† appears,
and we can calculate it using Eq. (<ref>).
Here, we use the following equations,

    ( 1 - γ_i )^†   = 1 - γ_i  ( i = 1,2,3 ) 
       = γ_5 ( 1 + γ_i ) γ_5  
    e^μ( 1 - γ_4 )^†   = e^μ^⋆( 1 - γ_4 )  
       = e^μ^⋆γ_5 ( 1 + γ_4 ) γ_5    .

Note that we now use the Euclidean gamma matrix {γ_i, γ_j } = 2 δ_i j.
Thus, we find

    [ Δ(μ) ]^⋆   = ( γ_5 Δ(-μ^⋆) γ_5 )  
       = ( γ_5 )  ( Δ(-μ^⋆) )  ( γ_5 )  
       = Δ(-μ^⋆)    .

Δ(μ) 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 [ Δ(μ) ]^N_f e^-S_G is complex.
Such a problem is called the sign problem<cit.>.
This problem was found in 1985<cit.>.
Although researchers have studied it about 30 years, it is not solved yet.









 §.§ Reweighting technique




The reweighting technique<cit.> may be the most famous method
which is made in order to avoid the sign problem.
In this method, we rewrite fermion matrix as

    Δ(μ) = Δ(μ)/Δ(μ_0)Δ(μ_0)    .

We usually set μ_0 = 0.
Then, partition function becomes

    Z(μ)    = ∫ DU [ Δ(μ)/Δ(μ_0)]^N_f[ Δ(μ_0) ]^N_f e^-S_G
       = Z(μ_0) ·1/Z(μ_0)∫ DU [ Δ(μ)/Δ(μ_0)]^N_f[ Δ(μ_0) ]^N_f e^-S_G
       = Z(μ_0) ·⟨[ Δ(μ)/Δ(μ_0)]^N_f⟩_μ_0  ,

where ⟨·⟩_μ_0 means the expectation value
on the gauge configurations which are generated at μ_0.
Similarly, the observables is written as

    ⟨ O(μ) ⟩_μ   = 1/Z(μ)∫ DU O(μ) [ Δ(μ)/Δ(μ_0)]^N_f[ Δ(μ_0) ]^N_f e^-S_G
       = 1/Z(μ) Z(μ_0) ·1/Z(μ_0)∫ DU O(μ) [ Δ(μ)/Δ(μ_0)]^N_f[ Δ(μ_0) ]^N_f e^-S_G
       = . ⟨ O(μ) [ Δ(μ)/Δ(μ_0)]^N_f⟩_μ_0/ ⟨[ Δ(μ)/Δ(μ_0)]^N_f⟩_μ_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 Carlo integration
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 Δ(μ) 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 effect of overlap problem.
In this method, we reweight e^-S_G too.
Concretely, 

    [ Δ(μ) ]^N_f e^-S_G(β) = [ Δ(μ)/Δ(μ_0)]^N_fe^-S_G(β)/e^-S_G(β_0)·[ Δ(μ_0) ]^N_f e^-S_G(β_0)  ,

where β is the effective coupling constant (see Eq. (<ref>), the definition of S_G).
Using this technique, we can enhance accuracy in part at high temperature<cit.>, T > T_c,
but the overlap problem still remains in the other part<cit.>.






 §.§ Taylor expansion method




We calculate the thermodynamic observables via the partition function Z(μ).
For example, the pressure is given as

    Δ p(μ, T)/T^4   := p(μ, T) - p(0, T)/T^4
       = 1/V T^3[ log(Z(μ, T)) - log(Z(0, T)) ]  
       = N_t^3/N_x N_y N_zlogZ(μ, T)/Z(0, T)  .

where p is the pressure and V = N_x a × N_y a × N_z a is the spatial volume of
the system.
Then, we can expand it in powers of μ/T around μ = 0,

    Δ p(μ, T)/T^4 = ∑_n=2,4,6,⋯^∞ c_n(T) (μ/T)^n    .

Here, the even powers of μ/T only survive because the QCD partition function has the

symmetry Z(μ) = Z(-μ).
We can calculate c_n via Z(μ) as follows,

    c_n    = 1/n !. ∂^n/∂ (μ/T)^nΔ p(μ, T)/T^4|_μ=0
       = 1/n !N_t^3/N_x N_y N_z. ∂^n/∂ (μ/T)^nlog Z(μ, T) |_μ=0  .

Because c_n does not depend on μ, we can calculate Δ p /T^4 at any μ using Eq. (<ref>)
once we calculate c_n.

In order to use such a method,
we have to calculate Z(μ) at several values of μ.
Although it is difficult to perform the calculation at finite μ because of the sign problem,
it is easy to perform near μ=0.
Thus, we calculate c_n 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 not work well at μ/T > 1 because the numerical errors become severe.
In fact, the studies of Taylor expansion method<cit.>
show the results of equation of state for μ/T ≤ 1.
(Several studies<cit.> show it for μ/T ≤ 1.2.)


Here, K. Nagata and A. Nakamura note that Taylor expansion and MPR methods
work out at the same region of μ<cit.>.
They concluded that they are consistent methods and the overlap problem of MPR
and the truncation 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.







 §.§ Canonical approach






Canonical approach is one method to avoid the sign problem.
It was suggested by A. Hasenfratz and D. Toussaint<cit.> in 1992.

They focused on the relation between grand partition function Z(μ)
and canonical partition function Z_C(n).
In the statistical mechanics,
it is well known that we can write Z(μ) as a polynomial of fugacity,

    Z(μ)    = Tr( e^-(Ĥ - μN̂) / T)  
       = ∑_n=-∞^∞⟨ n | e^-Ĥ/T| n ⟩ e^n μ/T
       =: ∑_n=-∞^∞ Z_C(n) (e^μ/T)^n    ,

with a thermodynamical limit.
Here, |n⟩ is the eigenvectors of the number operator N̂
(N̂ |n⟩ = n |n⟩) and Z_C(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 (<ref>) is called the fugacity expansion.

Hasenfratz and Toussaint noticed that Z(μ) is just the Laplace transform of Z_C(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 Z_C(n),

    Z(μ = iμ_I) = ∑_n=-∞^∞ Z_C(n) e^i (μ_I/T) n  .

Therefore, we also get Z_C(n) as the Fourier transform of Z(μ),

    Z_C(n) = 1/2 π∫ d( μ_I/T) Z(μ=iμ_I) e^-i(μ_I/T)n  .


Note that, for pure imaginary chemical potential,
Δ is real and we can calculate Z(iμ_I) and Z_C(n) without the sign problem.
Moreover, if we get Z_C(n), we can calculate Z(μ) with ANY chemical potential using Eq. (<ref>).
Thus, adopting such a procedure, we can calculate Z(μ) without the sign problem.
This method is called the canonical approach.

In this thesis, we use the canonical approach in order to avoid the sign problem.








§ REVIEW FOR STUDY OF QCD PHASE STRUCTURE WITH CANONICAL APPROACH






 §.§ Experiment data and canonical partition function


The canonical partition function Z_C(n) is related with the number of events
of the collider experiments.
This is because Z_C(n) represents the existence probability of the system
which n particles are in.
That is, we can get the following equation from Eq. (<ref>),

    1 = ∑_nZ_C(n) (e^μ/T)^n/Z(μ)  .

From this equation, we can consider Z_C(n) (e^μ/T)^n / Z(μ)
as the existence probability.





Here, let us consider the ratio Z_C(n) / Z_C(0).
This ratio is the normalized canonical partition function.
We will use it when we discuss the existence probability because

    . Z_C(n) (e^μ/T)^n/Z(μ)/ Z_C(0) (e^μ/T)^0/Z(μ). = Z_C(n)/Z_C(0)(e^μ/T)^n    .

Thus, we can cancel Z(μ) using the ratio.



We can see normalized existence probability Z_C(n) (e^μ/T)^n / Z_C(0) in collider experiments.
For example, Fig. <ref> shows the net proton dependence on the number of events
in the Au+Au collision at RHIC.


Let N_event(n_p) is the number of event with the net proton number n_p = 3 n
and N_total is the total number of events.
Then, normalized number of events N_event(n_p) / N_event(0) is given as follows,

    N_event(n_p)/N_event(0)   =  N_total·.Z_C(n_p) (e^μ/T)^n_p/ Z(μ). / N_total·.Z_C(0) (e^μ/T)^0 / Z(μ). 
       = Z_C(n_p)/Z_C(0)(e^μ/T)^n_p  .

Note that Z_C(n) satisfies

    Z_C(n) = Z_C(-n)    .

This is because the QCD partition function has the symmetry of the exchange
of perticles for anti–particles (C–symmetry).


Using Eqs. (<ref>) and (<ref>),
we can calculate Z_C(n) / Z_C(0) from RHIC data.
Figure <ref> shows the net proton number dependence on 
normalized canonical partition function from RHIC data.

From this figure, we see that Z_C(n) becomes smaller rapidly with increasing |n|
and it is raised with increasing the collision energy √(s_NN).







 §.§ 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<cit.>,
Hasenfratz and Toussaint performed the calculation of number density
on 2^4 and 4^4 lattice. See Fig. <ref>.
Their results tell us that, at low temperature, we can not perform this calculation well,
although at high temperature, we may perform it well.




In 2004, Y. Sasai, A. Nakamura and T. Takaishi investigated
the crossover phase transition of QCD at finite density<cit.>.
In this study, they used 8^3 × 4 lattice.
They calculated the chiral condensate as Fig. <ref>.

From this figure, we see the crossover phase transition as the jump of
the chiral condensate.
Thus, we can conclude that the canonical approach can be used in the study of
phase transition.
Note that, if T_c ≈ 200 MeV, they studied at μ≈ 100 MeV.


Ph. de Forcrand and S. Kratochvila studied the QCD phase diagram
via baryon number density 
with canonical approach in 2005<cit.>.
They then used 6^3 × 4 lattice.
Their results are shown in Figs. <ref> and <ref>.


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 are affected by the finite size effect strongly.
Therefore, we have to use larger lattice in order to investigate the QCD phase transition.




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. <cit.> in 2008. This is introduced below in detail.
Using this method, we can perform the calculation of the fermion determinant faster
than the other methods.



J. Danzer and C. Gattringer performed the calculation of number density
and number susceptibility at finite chemical potential with 8^3 × 4 lattice
in 2011<cit.>.
Then, they compare the results of canonical approach with
the results of Taylor expansion method as Fig. <ref> and Tab. <ref>.


From Fig. <ref>, we see that the results of canonical approach and Taylor expansion
method are consistent.
However, from Tab. <ref>,
we see the difference 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 even at this volume.
Since the fugacity series and Taylor series have different properties,
they suggested that we have to choose a better method (fugacity expansion or
Taylor expansion) for our purposes.














§ CALCULATION 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 very large.
Concretely, rank Δ = N_x N_y N_z N_t × 4 × 3
(N_x is a number of points of x coodinate, and so on.)
and if we choose N_x = N_y = N_z = N_t = 10, rank Δ≈ 10^5.
Moreover, we have to calculate it for EACH configuration. 





 §.§ Reduction formula





We firstly introduce the reduction formula of Wilson fermion matrix.
It was suggested by Nagata and Nakamura<cit.> in 2010.
(See also Ref.<cit.>.)
Note that the reduction formula for staggered fermion have been suggested since 1986
<cit.>.
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 - 2 z^-1κ r_- V - 2 z κ r_+ V^†  .

Here, r_± = (1 ±γ_4)/2 is the projection operators, z = e^-μ a and

    B(x, x^')    := δ_x, x^'
    		- κ∑_i=1^3( ( 1 - γ_i ) U_i (x) δ_x^', x+î
    			+ ( 1 + γ_i ) U_i^† (x^') δ_x^', x-î)  
         - κ C_SWδ_x, x^'∑_μ < νσ_μν F_μν (x)    ,  
    
    	V(x, x^')    := U_4 (x) δ_x^', x+4̂  ,  
    
    	V^†(x, x^')    := U_4^† (x) δ_x^', x-4̂  .

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

    B = 
    	[                                                                                                                                                                                                                                                                                                                                                                                             t^'=1 ⋯ t^'=N_t;                                                                                                                                                                 [   t=1;   t=2;   t=3;     ⋮; t=N_t ] (
    		[   B_1     0     0     ⋯     0     0;     0   B_2     0     ⋯     0     0;     0     0   B_3     ⋯     0     0;     ⋮     ⋮     ⋮     ⋱     ⋮     ⋮;     0     0     0     ⋯     0 B_N_t ])    , ]

















    V = (
    		[            0     U_4(t=1)            0            ⋯            0            0;            0            0     U_4(t=2)            ⋯            0            0;            ⋮            ⋮            ⋮            ⋱            ⋮            ⋮;            0            0            ⋯            ⋯ U_4(t=N_t-2)            0;            0            0            ⋯            ⋯            0 U_4(t=N_t-1);   U_4(t=N_t)            0            ⋯            ⋯            0            0;              ])   .



In order to cancel V^†, we multiply Δ and

    P := r_- + z^-1 r_+ V    .

As a result, 

    Δ P = (B r_- - 2 κ r_+) + (B r_+ - 2 κ r_-) z^-1 V    ,

Here, (B r_- - 2 κ r_+) is the diagonal part
and (B r_+ - 2 κ r_-) z^-1 V is the other part.
As the time–plane block matrix, it is given as 

    Δ P =
    		(
    		[          α_1     β_1 z^-1             ;                       α_2     β_2 z^-1             ;                         ⋱            ⋱             ;                   α_N_t-1 β_N_t-1 z^-1;   β_N_t z^-1        α_N_t ])    ,

where,

    α_i    = α^a b, μν(x⃗,y⃗,t_i)  
       = B^a b, μσ(x⃗,y⃗,t_i) r_-^σν
    			- 2 κ r_+^μνδ^a bδ(x⃗-y⃗)    ,  
    β_i    = β^a b, μν(x⃗,y⃗,t_i)  
       = B^a c, μσ(x⃗,y⃗,t_i) r_+^σν U_4^c b(y⃗,t_i)
    			- 2 κ r_-^μνδ(x⃗-y⃗) U_4^c b(y⃗,t_i)    .

Here, a, b is color indices

Using these equations, we can write the determinant of Δ P,

    (Δ P) = ( ∏_i=1^N_t(α_i) ) ( 1 + z^-N_t Q )    ,

where Q := (α_1^-1β_1) ⋯ (α_N_t^-1β_N_t).
Since (P) = z^-N/2 and N = N_x N_y N_z N_t × 3 × 4, we get

    Δ   = z^N/2( ∏_i=1^N_t(α_i) ) ( 1 + z^-N_t Q )  
       = z^-N/2( ∏_i=1^N_t(α_i) ) ( z^N_t + Q )    .

Then, we use a property of determinant (A B) = (A) (B),
and (cA) = c^N(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 is written as

    (z^N_t + Q) = ∏_n=1^N_red (λ_n + z^N_t)    ,

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

    z^-N/2∏_n=1^N_red (λ_n + z^N_t)
    		= ∑_n=-N_red/2^N_red/2 c_n ( z^N_t)^n    .

Here, using z^N_t = e^μ a N_t = e^μ/T, we finally obtain the following relation,

    Δ(μ) = ∑_n=-N_red/2^N_red/2 C_n ( e^μ/T)^n    .

We can calculate C_n numerically.
Nagata and Nakamura noticed later that C_n is just the canonical partition function Z_C(n).
Thus, we can calculate Z_C(n) with the eigenvalues of 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Δ = N_x N_y N_z N_t × 4 × 3
to rankα_i = rankβ_i = N_x N_y N_z × 4 × 3 = (rankΔ)/N_t.
The time complexity of a determinant calculation depending on the rank N_rank
is O(N_rank^3) in a typical algorithm.
In the reduction formula, we repeat the calculation of α_i and β_i N_t times.
Therefore, the complexity of its algorithm is O ( (N_x N_y N_z × 4 × 3)^3 × N_t ).
On the other hand, the complexity of naive method is O ( (N_x N_y N_z N_t × 4 × 3)^3 ).
If we set N_x = N_y = N_z = N_t = 10, it is reduced from O(10^15) to O(10^13).
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 on N_x, N_y, N_z and N_t
because the physics of system depends on them.

Although the reduction formula also seems efficient, calculation cost is still large.
Especially, when we calculate with large lattice,
the time complexity becomes to diverge with O ( (N_x N_y N_z)^3 ).
Therefore, we need to reduce it.







 §.§ 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:

    Δ   = exp[ TrlogΔ]  
       = exp[ Trlog( 1 - κ Q ) ]   .

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

    Δ (μ) = exp[ Tr∑_k=1^∞1/k(κ Q(μ) )^k]   .

This is the hopping parameter expansion.
Using this equation,
we can calculate Δ(μ) by the polynomial of Q(μ).


Then, the complexity is O ( (N_x N_y N_z N_t)^2).
This is because the trace of Q^k is calculated by multiplication of Q and a vector.
Such computation costs us O(N_rank^2) operations.
We can reduce the time complexity using this method.

However, the calculation cost is still large.
It is because we have to calculate Δ(μ) for varying μ.
In order to reduce the cost, we can use the following method.






 §.§ Winding number expansion




Winding number expansion (WNE)<cit.> is an improved method of HPE.

In HPE, we rewrite Δ to exp[ TrlogΔ]
and expand logΔ.
Then, using a method of statistical mechanics, we can rewrite the trace as

    TrlogΔ(μ)
    		= ∑_Ψ⟨Ψ| logΔ(μ) | Ψ⟩  ,

where | Ψ⟩ is a state in the system.
The right hand side of this equation means a summation 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.

Here, there are many types of size or shape of loops.
We therefore have to distinguish them to calculate in practice.
Now, let us focus on the chemical potential dependence.

From the definition of Δ(μ),
we notice that the coefficient exp(μ a N_t) = exp(μ/T) appears
when a fermion hop N_t times along the time direction.
We hence can write Δ(μ) formally,

    Δ(μ) = exp[ ∑_k=-∞^∞ W_k ( e^μ/T)^k ]   .

We can calculate W_k numerically. The algorithm is discussed in Appxs. <ref> and <ref>.
(See also Refs. <cit.> and <cit.>.)
Using this equation, we can see the μ dependence on Δ clearly.

Note that Eq. (<ref>) can be used at any chemical potential.
That is, we can calculate Z(μ) from Δ(μ) via WNE without canonical ensembles.
Such a calculation is discussed in Appx. <ref>.


Since W_k does not depend on μ, if ONCE we measure W_k,
we can calculate Δ with arbitrary μ.
The time complexity is

O ( (N_x N_y N_z)^2 ×. . (N_t)^2 )
in this method.
(See Appxs. <ref> and <ref>.)
The complexity of spatial part of this method is smaller than the one of reduction formula.
Therefore, we can calculate Z_C(n) in large lattice using WNE.
In this thesis, we adopt this method when we calculate Z_C(n).










§ PROBLEMS 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 one configuration Z^(1)_C(n).
We then generate a configuration at T/T_c = 0.93(5) and m_π/m_ρ = 0.80,
and set the sample size of Fourier transform as N_FT = 1024.

We find that Z_C(n) is a monotonic decrease function when n is small
and it becomes smaller rapidly with increasing n.
However, 

	
  * Z_C(n) is saturated at n = 50, and
	
  * it recovers near n = 1023 = N_FT-1.

As we know, Z_C(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.





 §.§ Property of discrete Fourier transform





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

DFT is a transformation which is defined as the following equation:

    F̃(k) = 1/N_FT∑_j=0^N_FT-1 f(j) e^-i ϕ_k j  ,
    		ϕ_k = 2 π k/N_FT  .

Now, we only have f(j) (j ∈𝐙 , 0 ≤ j < N_FT)
which are sampled on N_FT 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 + N_FT  n] j )    ,

with any n ∈𝐙.
That is why we can not distinguish F̃(k) from F̃(k+N_FT  n).
This phenomenon 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̃(N_FT) 
    F̃(-1)    = F̃(N_FT-1) 
    F̃(-2)    = F̃(N_FT-2)   .

Thus, Fig. <ref> shows that
n=1 system and “n=-1" system exists with the same probability when μ = 0.
This means that the symmetry of the exchange of particles for anti–particles
(C–symmetry) is satisfied in the QCD system.
Then, the highest baryon number is n_max = N_FT/2.

Such a property is seen in the work of H. Nyquist<cit.> in 1928
and proved by C.E. Shannon<cit.> in 1949.
n_max = N_FT/2 is called Nyquist frequency
or the folding frequency.


We may see this property from Fig. <ref>.
We noticed that Z_C(n) in Fig. <ref> is symmetric around n = N_FT/2 = 512.
Figure <ref> shows it clearly.

Note that this data is not symmetric completely
because there is the numerical errors which we will discuss below in detail.

Therefore, although we calculate Z_C(n) using DFT with N_FT, 
we can use only half of data.






 §.§ ^*An instability of DFT






We do not understand yet why Z_C(n) is saturated at n = 50.
Although it had not been investigated for a long time,
the author studied it from a numerical standpoint recently<cit.>. 


First, the author considered that this saturation is caused by the numerical
errors, and if we understand how the errors affect to our calculation, we can reduce this saturation.
He notes that we can classify the numerical errors into four types of errors<cit.>
(see also Refs.<cit.> and <cit.>).

	
  * 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 difference of mathematical and numerical values is called rounding error.

		For example, this error arises in the following type of calculation,
		
    1/3   = 3.333333⋯× 10^-1(Exact mathematical value)
        3.333333 × 10^-1(Numerical value; 7 digits.)  ,

		In this case, the exact mathematical value has infinite significant digits, 
		but the numerical value has only seven significant digits.
		The numerical result is less than the exact mathematical one as 3.333⋯× 10^-8.
		To perform rounding exactly, computers usually have a guard digit<cit.>.
			
	
  * Truncation error:
		when we estimate an infinite series numerically,
		we have to truncate it because we can use only finite computation resource.
		Then, there is a difference between the exact mathematical result and its truncated one.
		This difference is called truncation error.
		We can estimate the effect of this error by varying a truncation point of the series.


	
  * Cancellation of significant digits:
		when we calculate a subtraction between two nearly equal values,
		a lot of higher significant digits are canceled.
		Therefore, the result has only a few significant digits
		because the number of significant digits is limited in the numerical calculation.
		This phenomenon is called the cancellation of significant digits.

		It occurs in the following type of calculation:
		
    1.234567 - 1.234566 (7 digits)
        = 0.000001 (1 digit.)  .

		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 effect 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.)

		Then, although six significant digits are similarly lost in this calculation,
		sixteen digits still remain in the final result.


	
  * Loss of trailing digits:
		when we calculate an addition or subtraction between a huge number and small one,
		many lower significant digits of the small number are cut off.
		This is because the variables hold only finite significant digits.
		Then, the numerical result is underestimated. 
		This phenomenon is called the loss of trailing digits.

		Specifically, it arises the following type of calculation:
		
    3 × 10^10 + 2 × 10^0
        = 3.000000002 × 10^10(Exact mathematical result)
        3.000000 × 10^10(Numerical result; 7 digits.)   ,

		In this case, variables have seven significant digits,
		and the numerical result is underestimated from the mathematical one by 2 × 10^0.

		As the cancellation of significant digits, 
		we should increase the number of significant digits in order to reduce this loss.
		
    3 × 10^10 + 2 × 10^0
        = 3.000000002 × 10^10(Exact mathematical result)
        3.0000002000000 × 10^10(Numerical result; 16 digits.)   ,

		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
		since computers have a guard digit<cit.>.
		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 the transmission of rounding error as an effect of
		the cancellation or loss.

	
  * The truncation error may occur in our program.
		When we use HPE or WNE in order to calculate Δ(μ),
		we have to truncate the infinite series of these expansion.
		Although the truncation error arises in our calculation,
		we can control it and estimate its effect 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 found that the cancellation of significant digits occurs and cause the instability of DFT.
This cancellation arises almost at the last addition of the summation in Eq. (<ref>);
specifically,

    Z_C(n)    = 1/N[ ( ∑_k=0^N-2 Z ( 2 π k/N) e^-i (2 π k / N) n) .  
        + Z ( 2 π (N-1)/N) e^-i (2 π (N-1) / N) n] ∑_k=0^N-2.    ,  
    
       ↑This addition.


Figure <ref> shows how many digits are lost at the summation of Eq. (<ref>).
(Note that the author then uses the multi precision calculation which is described in section below.)
The algorithm is discussed in detail in Appx. <ref>.


From this figure, we can see that this cancellation becomes serious with increasing n.
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 the physical quantities,
but it has only 16 significant digits.
Thus, we can not calculate Z_C(n) well using the usual (double precision) variables.





 §.§ ^*Reduce these problems




In order to reduce these problems,
the author suggests to calculate Z_C(n) with large enough data space,
specifically,

	
  * N_FT 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 in the discussion below,
		whose size is 8^3 × 4.
		One considers that QCD critical point may be at roughly T_c ≈ 160 MeV and
		μ^c = μ^c_B/3 ≈ 200–400 MeV<cit.>.
		Therefore, at T = T_c, lattice spacing a is about 0.3 fm
		because N_t 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 N_FT/2 > 20.

	
  * 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 N_FT is large, the cancellation is serious,
		and then we have to increase 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 to estimate 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 Z_C(n) using the multi precision calculation.
This calculation extends the number of significant digits.
We can make it by binding several variables.
For instance, a real number x represents as follows in the base–b numerical system, 

    x = x_sign·( x_1 b^-1 + x_2 b^-2 + ⋯ + x_f b^-f) · b^x_0  ,

where

	
  * x_sign is the signature of x ( x_sign∈{ +1, -1 });

	
  * x_0 is the exponent of x (x_0 ∈𝐙);
	
	
  * x_1, x_2, ⋯, x_f are the mantissas (or significands,
			x_1 ∈{ 1, 2, ⋯, b-1 },
			x_i ∈{ 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, x_i are just the bits, x_0 is a set of several bits.
However, if we program this framework as a software (x_i are constructed 
by integer variables with any f, also b is a variable), we can change b and f,
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 Z_C(n) changing the precision of variables, 
and compared them as Fig. <ref> shows.












We find that, by increasing the number of significant digits f,
we can calculate Z_C(n) up to large n.
This is because Z_C(n) is just a remainder of the calculation of Eq. (<ref>).
That is, the number of canceled significant digits is directly related with a magnitude of Z_C(n).



He notes that, when we change the lattice volume V or temperature T
fixing the number of significant digits, a saturation point (kink) of Z_C(n) moves.

Z_C(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 Z_C(n) decreases slowly with increasing n.
We can see this actually in Fig. <ref>.

In the other situation, if only V becomes large,
also Z_C(n) decreases slowly too
since the net quark number density n/V becomes small.
If volume is changed as V → V^' = 2V,
Z_C(n) corresponds to Z^'_C(n^' = 2n)
because they represent the same density system.
Thus, we can replace Z_C(1) → Z^'_C(2),
Z_C(2) → Z^'_C(4), ⋯,
therefore Z_C(n) decreases slowly.














§ STUDY OF FINITE DENSITY PHASE TRANSITION WITH CANONICAL APPROACH



As seen above, by introducing these methods,
we can avoid the sign problem and calculate Z_C(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.




 §.§ Thermodynamical observables and partition function


From grand partition function Z(μ),
we can calculate thermodynamical observables as follows,

    p(μ_B,T)/T^4   = 1/V T^3log Z(μ_B,T)  
       = N_t^3/N_x N_y N_zlog Z(μ_B,T)   ,

where p is the pressure, V = N_x a × N_y a × N_z 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,
N_t a = 1/T, with lattice spacing a.

Taking a deviation of the pressure, we can get a useful relation,

    Δ p(μ_B,T)/T^4   = p(μ_B,T)/T^4 - p(0,T)/T^4
       = N_t^3/N_x N_y N_zlog( Z(μ_B,T)/Z(0,T))    .

It is useful because the partition function is normalized by Z(0, T)
as the expectation values of observables.

Other observables, baryon number density n_B and baryon
number susceptibility χ are given by the differential of pressure,

    n_B(μ_B,T)/T^3   = ∂/∂ (μ_B/T)p(μ_B,T)/T^4  , 
    χ(μ_B,T)/T^2   = ∂^2/∂ (μ_B/T)^2p(μ_B,T)/T^4  .




Note that, the partition function has the following periodicity as Eq. (<ref>),

    Z( μ/T) = Z( μ/T + i2 π k/3)    ,

which is caused by RW symmetry.
Using this equation, we can rewrite the partition function as

    Z( μ/T) = 1/3[ ∑_k=0^2 Z( μ/T + i2 π k/3) ]    .

Then, the canonical partition function Z_C(n) becomes

    Z_C(n, T)    = 1/2 π∫_0^2 π d( μ_I/T) [ 1/3∑_k=0^2 Z( i μ_I/T + i2 π k/3) ] e^i (μ_I/T) n
       = 1/2 π∫_0^2 π d( μ_I/T) [ 1 + e^2 π i / 3 + e^4 π i / 3/3 Z( i μ_I/T) ] e^i (μ_I/T) n  .

Therefore, Z_C has the following property,

    Z_C( n ≠ 3 k ) = 0    ,

for any k ∈𝐙.
This property is called the triality.

Thus, the grand partition function is written as

    Z(μ_B) = ∑_B=-∞^∞ Z_C(B) ( e^μ_B/T)^B  ,

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 the results of MPR in order to discuss the validity range of canonical approach.
We use the results of MPR in Ref. <cit.>,
and we adopt the same numerical setup as Ref. <cit.>
when we calculate our observables of canonical approach.













In Ref.<cit.>, 
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 be under control.
Therefore, our work enables us to check consistency among our canonical approach, MPR and Taylor expansion method.
















 §.§ Numerical setup


We adopted the clover improved Wilson fermion action with N_f = 2 and
C_SW = (1 - 0.8412/β)^-3/4 evaluated by a one–loop perturbation theory<cit.>,
and Iwasaki gauge action<cit.>.
All simulations were performed on a
N_x × N_y × N_z × N_t = 8 × 8 × 8 × 4 lattice.
Number of flavor is set to N_f = 2.
We considered values of β = 2.00, 1.95, 1.90, 1.85, 1.80 and 1.70;
they correspond to T/T_c = 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 for each value of β
by following the line of constant physics for the case of m_π / m_ρ = 0.80,
as in Ref. <cit.>.
In Ref.<cit.>, the lined of constant T/T_c is determined by the relation,
T_c/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. <cit.>
to evaluate T_c/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 so that the simulation time was N_τ× dτ = 1.
We neglected the first 2000 trajectories for termalization,
after we extracted 400 configurations every 200 trajectory for each parameter set.

Number of Fourier series N_FT is set to N_FT = 1024 for test in this thesis.
When we calculate the grand partition function Z(μ) from the canonical partition function Z_C(n),
we have to truncate the series of fugacity expansion.
The truncation point B_max is listed in Tab. <ref>,
which is discussed below in detail.
Number of significant digits f is set to f = 400 in this thesis.








 §.§ Error estimation




When we calculate the fugacity expansion of the grand partition function Z(μ),
we must truncate the series as

    Z(μ_B) = ∑_B=-B_max^B_max Z_C(B) ( e^μ_B/T)^B  .

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 B_max
and check the convergence of results.
If the results are stable when B_max is varied,
we can trust them because they should not depend on B_max physically.

For a convergence check, we evaluate the relative error R as follows<cit.>,

    R(μ_B) := ⟨ O(μ_B) ⟩_B_max - ⟨ O(μ_B) ⟩_B_max-1/⟨ O(μ_B) ⟩_B_max  ,

where ⟨ O(μ_B) ⟩_B_max is the expectation value
with the truncation at B_max.

In this thesis, we require that the results have at least two significant digits.
Therefore, we trust the results only for

    R(μ_B) < 10^-3  .







 §.§ ^*Numerical results







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




  §.§.§ Canonical partition function and analysis of truncation point



Before we introduce our results of thermodynamic observables,
we show the results of canonical partition function.

Figure <ref> shows the baryon number dependence on the canonical partition function.
From this figure, we can see that Z_C(B) becomes smaller rapidly with increasing B
at both of temperatures below and above T_c.
Also, we see that Z_C(B) depends on temperature; it is raised with increasing temperature.
We can consider that they are because Z_C(B) represents the existence probability of the system
which B particles are in[
	In detail, Z_C(B) (e^μ/T)^B / Z(μ) roles the weight function
	as we saw in Sec. <ref>,
	and then it represents the existence probability of the system.
	When μ=0, Z_C(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.
].
That is, B=1 system is more likely than B=2 or another system when μ=0,
and B=2 system at high temperature is easier to arise
than such systems at low temperature.



In Fig. <ref>, the canonical partition function is cut at B_max.
We define the value of B_max using the following way.
Because of numerical and statistical errors,
sometimes Z_C(B) becomes negative, although it has to satisfy Z_C(B) > 0.
Then, if the error of Z_C(B) is larger than the absolute value of Z_C(B),
we can use the positive value -Z_C(B) instead of the negative value Z_C(B).
However, if Z_C(B) is negative and δ[ Z_C(B) ] < | Z_C(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 B_max from the previous data.







  §.§.§ Pressure, baryon number density and baryon number susceptibility


Let us introduce our results of thermodynamic observables.

First, we examine the pressure.
Fig. <ref> shows the chemical potential dependence on pressure.
From this figure, it is found that the results of pressure at high temperature T > T_c 
do not suffer from large errors up to μ_B/T = 5, approximately. 
The results at low temperature T < T_c 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 suffered from fluctuations caused by the (crossover) phase transition at zero density. 
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. <ref> 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 the results
of MPR method.







Next, we consider the expectation value of the baryon number density.
Fig. <ref> shows the chemical potential dependence on baryon number density.
This figure demonstrates that,
for temperatures above and below T_c, the results are reliable up to μ_B/T = 4
and μ_B/T = 3–3.5, respectively.
However, near T_c, the reliable chemical potential range are limited up to μ_B/T = 2.4.
This may be for the same reason described in the pressure analysis.

Figure <ref> represents the comparison of baryon number density between the methods.
It demonstrates good agreement between the results of the canonical approach
and those of MPR method in this case.

Moreover, the gradient of the baryon number density n_B/T^3 
as a function of baryon chemical potential μ_B
becomes small as the temperature decreases.
In a zero temperature case, 
n_B is expected to be zero up to μ_B=m_B where m_B is the lightest baryon mass of a system,
and becomes a finite value at this point.
This phenomenon is called the silver blaze phenomenon<cit.>.
The data at T/T_c=0.84 in Fig. <ref> does in fact show such a feature.






Finally, we investigate the baryon susceptibility.
Figure <ref> shows the results of the susceptibility.
From this figure, we find that the results are reliable up to μ_B/T = 3.5 at T > T_c,
whereas they are reliable up to T_c = 2.4–2.9 at T < T_c. 

Figure <ref> represents the comparison between the methods in the susceptibility case.
We see that the susceptibility results of canonical approach are in very good agreement
with those of the MPR method.

The baryon susceptibility as a function of μ_B/T does not show a clear peak.
Thus, signals of the finite density transition between the confinement–deconfinement phases
is not observed yet in this case.








§ CONCLUSION 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 of
finite density QCD because other methods, such as MPR, Taylor expansion and 
imaginary chemical potential method<cit.>, yield reliable information 
in practical situations only up to μ_B/T=3.
Here, the multi precision calculation significantly contributes to this conclusion.

We have investigated the QCD phase structure via the pressure, baryon number density
and baryon number susceptibility.
We calculated these thermodynamic observables in 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 more accurately at large baryon numbers.
As shown in Fig. <ref>,
our canonical partition functions currently have unphysical phases at several baryon numbers
although these should have a real positive value in principle.
After the author works, ones have studied
why such phases appear in the canonical partition function Z_C(n)<cit.>.
However, we have not understood it well.

In our opinion, this may be originated from the overlap problem<cit.>.
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 with gauge configurations generated 
at each pure imaginary chemical potentials to realize the appropriate 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 have been applied to the other works.
For example, the author and his collaborators calculated the higher cumulants
with our improvement of canonical approach<cit.>.
However, the error 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
Z_C(n) with more high accuracy<cit.>.
Then, they used fitting procedure for the number density at imaginary chemical 
potential and constructed Z_C(n) from it with more high accuracy.

In addition, they calculated the higher cumulants
with their procedure based on our improvement of canonical approach<cit.>.
Then, they showed that their results was consistent to experimental results.



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 of gauge configurations.

The author considers that
such methods may include the calculation which causes the cancellation 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 is always effectual for such methods.

The author notices that
we can calculate Z_C(n) with high enough accuracy using the reduction formula
even though we do not use the multi precision calculation[
	Then, we only need to compute the mantissas and exponent part respectively.
	This is because very small eigenvalues appear (λ < 10^-308) in this calculation.
	We do not need to extend the number of significant digits.
][
	But, as seen above, the calculation cost of reduction formula is much larger
	than the cost of winding number expansion and Fourier transformation.
].
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
<cit.>. 
We can also find through our work that the canonical approach 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. (<ref>) with high accuracy.
This requires more computational resources than what has been reported here
but is within the scope of the next–generation high–performance era. 


















§ 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 me beer, sake, wine, comics, PCs and money.



My calculations were performed
using the supercomputers, SX–9 and SX–ACE at Reserach Center for Nuclear Physics (Osaka Univ.),
PACS–IX and HA–PACS at Center for Computational Sciences (Univ. of Tsukuba),
SR16000 at Yukawa Institute for Theoretical Physics (Kyoto Univ.)
and FX10 at Information Technology Center (Univ. of Tokyo).

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].





















PART: 
*Appendix






§ GRASSMANN INTEGRAL OF FERMION ACTION




Fermion fields are the Grassmann numbers which satisfy the anticommutation relation,

    Ψ_i Ψ_j    = -Ψ_j Ψ_i    ,  
    Ψ̅_i Ψ̅_j    = -Ψ̅_j Ψ̅_i    .

The Grasmann integral is defined by the following relations,

    ∫ dΨ_i    = ∫ dΨ̅_i = 0    ,  
    ∫ dΨ_i  Ψ_i    = ∫ dΨ̅_i  Ψ̅_i = 1    .

Also, we get the following relations from Eq. (<ref>),

    Ψ_i Ψ_i = Ψ̅_i Ψ̅_i = 0    .

From this equation, we find that the polynomial of Grassmann numbers includes 
the terms 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(-Ψ̅𝐀Ψ) (𝐀 is a matrix)
is calculated as

    ∫ dΨ_1 dΨ̅_1 dΨ_2 dΨ̅_2   e^-∑_i,j=1^2Ψ̅_i A_i jΨ_j
       = ∫ dΨ_1 dΨ̅_1 dΨ_2 dΨ̅_2   e^-Ψ̅_1 A_1 1Ψ_1 e^-Ψ̅_1 A_1 2Ψ_2 e^-Ψ̅_2 A_2 1Ψ_1 e^-Ψ̅_2 A_2 2Ψ_2
       = ∫ dΨ_1 dΨ̅_1 dΨ_2 dΨ̅_2  ( 1 - Ψ̅_1 A_1 1Ψ_1 ) ( 1 - Ψ̅_1 A_1 2Ψ_2 )  
       = dΨ_1 dΨ̅_1 dΨ_2 dΨ̅_2  ×( 1 - Ψ̅_2 A_2 1Ψ_1 ) ( 1 - Ψ̅_2 A_2 2Ψ_2 )  
       = (-A_11) (-A_22) + (-A_21) (+A_12) = 𝐀  .

We used the Taylor expansion of the exponential in the third line.
Note that the higher powers (e.g. (Ψ̅_1 A_1 1Ψ_1)^2) do not appear
because of Eq. (<ref>).
We also use Eq. (<ref>) and Eq. (<ref>)
in the fourth line.

We can generalize this equation to more general formulation.
T. Matthews and A. Salam give us the following formulas<cit.>,

    ∫ DΨ DΨ̅  e^-Ψ̅ΔΨ   = Δ  ,  
    ∫ DΨ DΨ̅ Ψ_i Ψ̅_j e^-Ψ̅ΔΨ   = -( Δ^-1)_i jΔ  ,

where Δ is the fermion matrix.
These are called Matthews–Salam formula.







§ ^*ALGORITHM OF WINDING NUMBER EXPANSION








We can calculate the coefficient W_k of Eq. (<ref>) via HPE<cit.>.
This is because Tr Q^m is given as follows,

    Tr Q^m = ∑_x, a, α⟨Ψ^a, α(x) | Q^m | Ψ^a, α(x) ⟩  ,

where a and α denote the Dirac and color index, respectively.
Here, we dissolve Q into the spatial component Q_i (i=1, 2, 3),
forward temporal component Q_4^(+) and backward temporal component Q_4^(-),

    Q^m | Ψ^a, α(x) ⟩ = ( ∑_i=1^3 Q_i + e^+μ a Q_4^(+) + e^-μ a Q_4^(-))^m | Ψ^a, α(x) ⟩  .

Note that Q_i includes the clover term.
Then, we can classify this with chemical potential dependence.
That is,

    Q^m | Ψ^a, α(x) ⟩   = ∑_k=-m^+m( e^+μ a)^k Q^(k)(m) | Ψ^a, α(x) ⟩
       =: ∑_k=-m^+m( e^+μ a)^k X⃗^(k)(m; x, a, α)    ,

where Q^(k)(m) is the product of Q_i and Q_4^±.

Thus, we can write Tr Q^m as

    Tr Q^m = ∑_x, a, α⟨Ψ⃗(x, a, α), ∑_k=-m^+m( e^+μ a)^k X⃗^(k)(m; x, a, α) ⟩  ,

where Ψ⃗ := | Ψ⟩ and
⟨· , ·⟩ is the inner product.

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 ( (N_x N_y N_z N_t × 3 × 4)^2 )
to O ( (N_x N_y N_z × 3 × 4)^2 ).
In the practical computing, we do not perform the summation of m.
We thus save ⟨Ψ⃗, (e^μ a)^k X⃗^(k)⟩ =: Z^(k)
as an array for k.
After we perform the summation of m,
we extract the coefficient of (e^μ/T)^k = (e^μ a)^k N_t.
The coefficient is just W_k.



Since the algorithm is intricate, we show it as a pseudo code herein.
Note that the outside loops of (x,y,z,t,c,μ) for the trace
are performed by the noise method which is discussed below.








§ NOISE METHOD




The trace of matrix is performed by the loops of summation.
Naively, we perform the loops N_x N_y N_z N_t × 3 × 4 times.
We can reduce this calculation using the noise method<cit.>.

In this method, we start from Eq. (<ref>), that is,

    Tr𝐀 = ∑_ϕ⟨ϕ| 𝐀| ϕ⟩  ,

where 𝐀 is a matrix with the rank of N_rank.
Then, | ϕ⟩ is a vector.
Here, we use the notation ϕ⃗ := | ϕ⟩ as above section;
we then write as

    Tr𝐀 = ∑_ϕ⟨ϕ⃗, 𝐀ϕ⃗⟩  ,

where ⟨·, ·⟩ is the inner product.


We can calculate it choosing ϕ as the vector of noises<cit.>,

    Tr𝐀 = lim_N →∞1/N∑_i=1^N⟨η⃗^(i), 𝐀η⃗^(i)⟩  ,

where η⃗^(i) is the i–th noise vector.
The noise vector satisfies the following property,

    1/N_rank⟨η⃗^(i), η⃗^(j)⟩ = δ_i j  .

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
with high enough accuracy<cit.>.
In this thesis, for safety, we use 64 noise vectors for high temperature T > T_c simulations and
128 noise vector for low temperature T < T_c simulations.
An example calculation is shown as follows. See Fig. <ref>.

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 is because the times of loops are reduced as
N_x N_y N_z N_t × 3 × 4 = 24576 → N_noise = 128
when the lattice size is 8^3 × 4.
Then, the time complexity of WNE with noise trace is

O ( (N_x N_y N_z)^2 × (N_t)^2 ).










§ ^*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,

    W_k → Z(i μ_I) → Z_C(n) → Z(μ)    .

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 W_k directly,

    W_k → Z(μ)    ,

because we can use the WNE of fermion determinant at any μ (not only i μ_I).
(See Fig. <ref>.)

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 differs 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,
			
    Δ(μ) = exp[ ∑_k W_k (e^μ/T)^k ]   ,

			before we calculate the integral of gauge field.
			On the other hand, in the canonical approach,
			we perform the finite density calculation
			Z(μ) = ∑_n Z_C(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 [ Δ(μ) ]^N_f e^S_G
using complex Δ(μ).


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 <ref> shows the imaginary part of pressure
which we calculate using the canonical approach and direct method.

From this figure, we see that the results of canonical approach and direct method
are consistent to zero.
However, the statistical error of canonical approach is smaller than
the one of direct method.
We can see this property from other observables.
For instance, Fig. <ref> and <ref>
show the imaginary part of baryon number density and baryon susceptibility.
We find that the error of canonical approach is always smaller than the one of direct method.















Figure <ref>, <ref> and <ref> show the thermodynamic observables
which we calculate using the canonical approach and direct method.
Then, the validity range of canonical approach is estimated as Sec. <ref>.







The range of direct method is estimated
by the magnitude of statistical error of imaginary part.
The procedure is as follows.




	
  * Save the error of imaginary part of observables of canonical approach
		at the boundary of validity range, e.g.
		δ( Im[ Δ p(μ_val(can)) / T^4 ]_can).
	
  * Compare the one of direct method
		δ( Im[ Δ p(μ) / T^4 ]_dir) with

		δ( Im[ Δ p(μ_val(can)) / T^4 ]_can).
	
  * Find μ_1 which satisfies,
		
    δ( Im[ Δ p(μ_1) / T^4 ]_dir)
    			= δ( Im[ Δ p(μ_val(can)) / T^4 ]_can)    .

	
  * 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 than the one of direct method.
Especially, in the pressure data, canonical approach is valid beyond μ_B/T=5 at T/T_c=1.20,
although direct method is valid to μ_B/T=4.

Therefore, when we calculate the thermodynamic observables at finite density,
it is better that we use the canonical approach in order to avoid the sign problem.







§ ^*ALGORITHM FOR ESTIMATION OF CANCELLATION OF SIGNIFICANT DIGITS




In 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, the logarithm of a = 6.022 × 10^23 is log_10 a = 23.78.
Then, “23" corresponds to the power of a,
and “0.78" corresponds to the significant digits of a.
This is because we can dissolve log_10 a as

    log_10 a = log_10 (6.022 × 10^23) = log_10 (6.022) + log_10 (10^23)    .

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    ,  
    
     	d    := log_10 (a) - log_10 (r)    .

where “r" means the result and “d" means the difference between the magnitudes of a and r.
If a and b are

    a    = 6.022 × 10^23  ,  
    
    	b    = 6.021 × 10^23  ,

then r and d become the following numbers,

    r    = log_10 (0.001 × 10^23) = 20.00    ,  
    
     	d    = 23.78 - 20.00 = 3.78    .


In this case, when we calculate r = a - b, three significant digits are lost.
This “three" appears in Eq. (<ref>).
That is, the number of canceled digits are given as

    ( Number of canceled digits) = ⌊ d ⌋ = ⌊log_10 (a) - log_10 (r) ⌋  ,

where ⌊·⌋ is the floor function.

Here, what does “0.78" of Eq. (<ref>) mean?
This is the difference between the magnitudes of significant digits.
Let a_sig and a_pow be the significant digits and power of a, respectively.
We then can write d as

    d    = log_10 (a) - log_10 (r) = log_10 (a_sig× 10^a_pow) - log_10 (r_sig× 10^r_pow)  
       = log_10( a_sig/r_sig) + log_10( 10^a_pow/10^r_pow)  
       = log_10( a_sig/r_sig) + ( a_pow-r_pow)    ,  
    ⌊ d ⌋   = a_pow-r_pow  .

Thus, “0.78" is log_10( a_sig/r_sig).
(In fact, we can compute it as log_10 (6.022/ 1.0) = 0.78.)
Note that 1 ≤ a_sig/r_sig < 10.
This is because a_sig and r_sig satisfy
1 ≤ a_sig < b and 1 ≤ r_sig < b where b is the base of logarithm.
(b=10; See also Sec. <ref>.)


The Algorithm for estimation of the number of canceled digits is as follows.
Note that we now adopt the multi precision calculation.

	
  * Preparation: adopt the multi precision calculation to the program.

	
  * Initialization:
		
			
  * Set the number of significant digits is large enough.
			
  * Choose a as the minuend and b as subtrahend.
				They are multi precision variables.
		

	
  * Estimation:
		
			
  * Save log_10 (a).
			
  * Calculate the subtraction r = a - b which may cause the cancellation of significant digits.
			
  * Calculate d = log_10 (a) - log_10 (r).
			
  * Print ⌊ d ⌋ which is the number of canceled digits.
		










99



	
	
	Fukushima
		K. Fukushima and T. Hatsuda, “The phase diagram of dense QCD,"
		 	Rept. Prog. Phys. 74 (2011) 014001.


	QM
		M. Buballa, and M. Oertel, “Color–flavor unlocking and phase diagram with
		self–consistently determined strange-quark masses,"
		Nucl. Phys. A 703 (2002) 770.


	

	beta_function
		C. Itzykson, and J. B. Zuber, “Quantum Field Theory," Dover Publications (2006).

	SLAC-exp
		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,
		“High-Energy Inelastic e–p Scattering at 6^∘ and 10^∘,"
		Phys. Rev. Lett. 23 (1969) 930.
	
	SLAC-th
		J.D. Bjorken and E.A. Paschos,
		“Inelastic Electron–Proton and γ–Proton Scattering and the Structure of the Nucleon,"
		Phys. Rev. 185 (1969) 1975.

	Gross
		D.J. Gross, F. Wilczek. “Ultraviolet behavior of non-abelian gauge theories,"
		Phys. Rev. Lett. 30 (1973) 1343. 
	
	Politzer
		H.D. Politzer. “Reliable perturbative results for strong interactions,"
		Phys. Rev. Lett. 30 (1973) 1346. 

	Creutz-1
		M. Creutz, “Monte Carlo study of quantized SU(2) gauge theory," Phys. Rev. D 21 (1980) 2308.

	Creutz-2
		M. Creutz, “Asymptotic-freedom scales," Phys. Rev. Lett. 45 (1980) 313.
		

	

	finite_temperature_QCD
		D.J. Gross, R.D. Pisarski and L. G. Yaffe,
		“QCD and instantons at finite temperature,"
		Rev. Mod. Phys. 53 (1981) 43.

	RW
		A. Roberge and N. Weiss. “Gauge theories with imaginary chemical potential and the phases of QCD,"
		Nucl. Phys. B 275 (1986) 734.

	Fodor-Katz
		Z. Fodor and S.D. Katz,
		“Critical point of QCD at finite T and , lattice results for physical quark masses,"
		JHEP 0404 (2004) 050.


	Taylor_critical
		S. Datta, R.V. Gava and, S. Gupta,
		“The QCD Critical Point : marching towards continuum,"
		Nucl. Phys. A 904 (2013) 883.

	models
		O. Scavenius, A. Mocsy, I.N. Mishustin and D.H. Rischke,
		“Chiral Phase Transition within Effective Models with Constituent Quarks,"
		Phys. Rev. C 64 (2001) 045202.




	

	kappa
		K. Jansen and R. Sommer,
		“O(a) improvement of lattice QCD with two flavors of Wilson quarks,"
		Nucl. Phys. B 530 (1998) 185.

	Iwasaki
		Y. Iwasaki, “Renormalization Group Analysis of Lattice Theories and
		Improved Lattice Action. II — four–dimensional non–abelian SU(N) gauge model —,"
		UTHEP-118 (1983), preprint.


	clover_fermions
		B. Sheikholeslami and R. Wohlert,
		“Improved continuum limit lattice action for QCD with Wilson fermions,"
		Nucl. Phys. B259 (1985) 572.





	sign1
	A. Nakamura, “Quarks and Gluons at Finite Temperature and Density,"
	Phys. Lett. 149B (1984) 391.

	sign2
	Ph. de Forcrand, “Simulating QCD at finite density," PoS (LATTICE2009) 010 (2010).


	sign_orig1
		J. Engels and H. Satz,
		“Deconfinement at finite baryon number density," Phys. Lett. B 159 (1985) 151.
	
	sign_orig2
		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.

	

	reweighting-1
		Z. Fodor and S.D. Katz, “A new method to study lattice QCD at finite temperature
		and chemical potential," Phys. Lett. B 534 (2002) 87.

	reweighting-2
		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.

	NN_eos
		K. Nagata and A. Nakamura, “EoS of finite density QCD with Wilson fermions
		by Multi–Parameter Reweighting and Taylor expansion," JHEP 04 (2012) 092.


	

	Taylor-1
		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.


	Taylor-2
		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 to 𝒪(μ_B^ 6) from Lattice QCD,"
		Phys. Rev. D 95 (2017) 054504.

	Taylor-3
		C. Schmidt and S. Sharma, “The phase structure of QCD," J. Phys. G 44 (2017) 104002.


	WHOT
		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.


	

	canonical
		A. Hasenfratz and D. Toussaint, “Canonical ensembles and nonzero density quantum chromodynamics,"
		Nucl. Phys. B 371 (1992) 539.




	RHIC
		M.M. Aggarwal, et al,
		“Higher Moments of Net Proton Multiplicity Distributions at RHIC,"
		Phys. Rev. Lett. 105 (2010) 022302.


	RHIC_fit
		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).


	SNT
		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


	Forcrand-1
		S. Kratochvila and Ph. de Forcrand,
		“The canonical approach to Finite Density QCD,"
		PoS (LATTICE2005) 167 (2006).
		
	Forcrand-2
		Ph. de Forcrand and S. Kratochvila,
		“Finite density QCD with a canonical approach,"
		Nucl. Phys. B (Proc. Suppl.) 153 (2006) 62.


	WNE
		X. Meng, A. Li, A. Alexandru and K.F. Liu,
		“Winding number expansion for the canonical approach to finite density simulations,"
		PoS (LATTICE2008) 032 (2009).


	Gatt
		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.




	NN-reduction
		K. Nagata and A. Nakamura, “Wilson Fermion Determinant in Lattice QCD,"
		Phys. Rev. D 82 (2010) 094027.

	reduction-2
		A. Alexandru and U. Wenger, “QCD at non–zero density and canonical partition
		functions with Wilson fermions,"
		Phys. Rev. D 83 (2011) 034502.
		
	reduction_stag
		P.E. Gibbs,
		“The fermion propagator matrix in lattice QCD,"
		Phys. Lett. B 172 (1986) 53.

  

	WNE_Gatt
		J. Danzer and C. Gattringer,
		“Winding expansion techniques for lattice QCD with chemical potential,"
		Phys. Rev. D 78 (2008) 114506.
  


	Oka-3
		A. Nakamura, S. Oka and Y. Taniguchi, 
		“QCD phase transition at real chemical potential with canonical approach,"
		JHEP 02 (2016) 054.


	Oka-4
		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).





	

	Nyquist
		H. Nyquist, “Certain topics in telegraph transmission theory,"
		Trans. AIEE. 47 (1928) 617.

	Shannon
		C. E. Shannon, “Communication in the Presence of Noise,"
		Proceedings of the IRE 37 (1949) 10.


	

	Oka-1
		S. Oka, “Exploring finite density QCD phase transition with	canonical approach
		—Power of multiple precision computation—,"
		PoS (LATTICE2015) 067 (2016).

	Oka-2
		R. Fukuda, A. Nakamura and S. Oka, “Canonical approach to finite density QCD
		with multiple precision computation,"
		Phys. Rev. D 93, 094508 (2016).


	TheArt
		D. E. Knuth,
		“The Art of Computer Programming, Volume 2 / Seminumerical Algorithms,"
		Second Edition, Addison–Wesley publishing (1969).

	guardbit
		D. Goldberg, “What every computer scientist should know about floating-point arithmetic,"
		ACM Computing Surveys (CSUR) 23 (1991) 5.












	
	
	WHOT_ref
		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.

	Oka-5
		R. Fukuda, A. Nakamura and S. Oka,
		“Validity range of canonical approach to finite density QCD,"
		PoS (LATTICE2015) 167 (2016).


	

	Gatt-SB
		C. Gattringer and T. Kloiber,
		“Lattice study of the Silver Blaze phenomenon for a charged scalar ϕ^4 field,"
		Nucl. Phys. B 869 (2013) 56.


	



	imag_chem
		Ph. de Forcrand and O. Philipsen,
		“The QCD phase diagram for small densities from imaginary chemical potential,"
		Nucl. Phys. B 642 (2011) 290.




	Suzuki-0
		R. Fukuda, A. Nakamura, S. Oka, A. Suzuki and Y. Taniguchi,
		“Study of the sign problem in canonical approach,"
		PoS (LATTICE2016) 060 (2017).

	Russia-0
		V.A. Goy, V. Bornyakov, D. Boyda, A. Molochkov, A. Nakamura, A. Nikolaev and V. Zakharov,
		“Sign problem in finite density lattice QCD,"
		PTEP 2017 (2017) 031D01.

	Suzuki-1
		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).

	Russia-1
		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.

	Russia-2
		V.G. Bornyakov, D.L. Boyda, V.A. Goy, A.V. Molochkov, A. Nakamura, A.A. Nikolaev and V.I. Zakharov,
		“New approach to canonical partition functions computation in N_f=2 
		lattice QCD at finite baryon density,"
		Phys. Rev. D 95 (2017) 094506.
	
	Russia-3
		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.




	Alex
		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.

	Li:2010qf 
		A. Li, A. Alexandru, K.F. Liu and X. Meng,
		“Finite density phase transition of QCD with N_f=4 and N_f=2 using canonical ensemble method,"
		Phys. Rev. D 82 (2010) 054502.

	Li:2011ee 
		A. Li, A. Alexandru and K.F. Liu,
		“Critical point of N_f = 3 QCD from lattice simulations in the canonical ensemble,"
		Phys. Rev. D 84 (2011) 071503.

	Gat1 
		A. Alexandru, C. Gattringer, H.P. Schadler, K. Splittorff and J.J. M. Verbaarschot, 
		“Distribution of canonical determinants in QCD,"
		Phys. Rev. D 91 (2015) 074501.
  
	Gat2 
		C. Gattringer, H.P. Schadler, 
		“Generalized quark number susceptibilities from fugacity expansion at finite chemical
		potential for N_f = 2 Wilson fermions,"
		Phys. Rev. D 91 (2015) 074511.



	Matthews-Salam
		T. Matthews and A. Salam, “Propagators of Quantized Field," Nuovo Cimento 2 (1955) 120.



	noise
		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.


















empty
top=2.7cm, bottom=0.5cm, left=0.5cm, right=0.5cm










Extra




In this section, there are some illustrations which are drawn by the author.






 	
    < g r a p h i c s >










































