Analytic Calculation of Ground State Properties of The 2D and 3D Electron Gas

ANALYTIC CALCULATION OF GROUND
STATE PROPERTIES OF THE 2D AND 3D

ELECTRON GAS
a thesis submitted to
the graduate school of engineering and science
of bilkent university
in partial fulfillment of the requirements for
the degree of
master of science
in
physics
By
Yağmur Katı
June, 2015
Analytic Calculation of Ground State Properties of the 2D and 3D
Electron Gas
By Yağmur Katı
June, 2015
We certify that we have read this thesis and that in our opinion it is fully adequate,
in scope and in quality, as a thesis for the degree of Master of Science.
Asst. Prof. Dr. Balázs Hetényi(Advisor)
Prof. Dr. Bilal Tanatar
Assoc. Prof. Dr. Hande Toffoli
Approved for the Graduate School of Engineering and Science:
Prof. Dr. Levent Onural

Director of the Graduate School
ii
ABSTRACT
ANALYTIC CALCULATION OF GROUND STATE
PROPERTIES OF THE 2D AND 3D ELECTRON GAS
Yağmur Katı
M.S. in Physics
Advisor: Asst. Prof. Dr. Balázs Hetényi
June, 2015
The electron gas (2D and 3D) is a model which consists of interacting electrons
moving in a uniform positive background. Its importance stems from the fact that
a number of metals behave similarly, it provides the functional used in density
functional theory, and that in 2D it can be experimentally realized. Understand-
ing the behavior of this model is of fundamental importance. In this thesis we
present an analysis of this model based on the Hypernetted Chain Method in 3D,
and 2D. The HNC method is a variational method to calculate the ground state
properties of an interacting system, by expressing the ground state energy as a
functional of the radial distribution function. Minimizing the energy expression
one obtains a zero energy Schrödinger equation for the square root of the radial
distribution function. The potential in this equation can include the effects of
fermionic or bosonic exchange. We applied this method to charged boson and
electron gas in 2D and 3D systems. On the basis of the results of this research, it
can be concluded that we obtained very close correlation energy results compared
to Monte Carlo, and FHNC results for the density range when rs is from 0 to 20.
This extended range is important for solid state applications.
Keywords: Electron gas, charged boson, quantum Hall effect, Hypernetted Chain
Theory, radial distribution function, ground state energy.
iii
ÖZET
İKİ VE ÜÇ BOYUTLU ELEKTRON GAZININ TEMEL
DURUM ÖZELLİKLERİNİN ANALİTİK
HESAPLAMASI
Yağmur Katı
Fizik, Yüksek Lisans
Tez Danışmanı: Asst. Prof. Dr. Balázs Hetényi
Haziran, 2015
İki ve üç boyutlu elektron gaz modeli, tekdüze pozitif yüklü bir alanda hareket
eden etkileşim halindeki elektronlardan oluşur. Birçok metal birbiriyle ben-
zer davranış gösterdiği, DFTde kullanılan fonksiyoneli sağladığı ve iki boyutta
hesaplanan sonuçları deneysel olarak da gerçekleştirilebildiği için bu modelin
davranışını anlamak temel önem taşımaktadır. Bu çalışma ile sizlere, HNC meto-
dunu kullanarak üç boyutlu ve iki boyutlu sistemlerde elektron gaz modelinin bir
analizini sunuyoruz. HNC Metodu, etkileşim halindeki bir sistemin temel durum
enerjisini radyal dağılım fonksiyonunun bir fonksiyoneli olarak ele alıp, temel du-
rum özelliklerini hesaplayan varyasyonel bir metoddur. Enerji ifadesi minimize
edilirse, radyal dağılım fonksiyonunun karekökü için sıfır enerjili Schrödinger den-
klemi elde edilir. Bu denklemdeki potansiyel, fermiyonik veya bozonik değişim
içerebilir. Biz bu yöntemi iki ve üç boyutlu yüklü bozon ve elektron gazına ve
quantum Hall sistemine uyguladık. Bu araştırmanın sonuçlarına dayanarak, katı
hal fiziği uygulamalarında önemli olan 0 < rs < 20 genişletilmiş yoğunluk aralığı
için Monte Carlo ve FHNC sonuçlarına çok yakın korelasyon enerjisi sonuçları
elde ettiğimizi söyleyebiliriz.
Anahtar sözcükler : Elektron gazı, yüklü bozon gazı, kuantum Hall etkisi, temel
durum enerjisi.
iv
Acknowledgement
The author wishes to thank several people. I would like to express my special
appreciation and thanks to my advisor Dr. Balázs Hetényi, for encouraging my
research and for allowing me to grow as a scientist.
I would also like to thank my committee members, Prof. Bilal Tanatar and
Dr. Hande Toffoli for their brilliant comments and suggestions; and I am grateful
to Prof. Cemal Yalabık, and Prof. Atilla Aydınlı for their generous guidance and
support.
I am hugely indebted and throughly grateful to Sinan Gündoğdu for his pre-
cious help and support during the time it has taken me to finalize this thesis.
Finally I would like to extend my deepest gratitude to my family for all of the
sacrifices that they have made on my behalf.
v
Contents
1 Introduction 1
2 Background 3
2.1 Jellium Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1.1 Particle Density and Wigner-Seitz Radius . . . . . . . . . 4
2.2 Coulomb Interaction . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Radial Distribution Function . . . . . . . . . . . . . . . . . . . . . 6
2.3.1 Static Structure Factor . . . . . . . . . . . . . . . . . . . . 8
2.4 Description of a Quantum Many Body System . . . . . . . . . . . 9
3 Hypernetted Chain Method 13
3.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Minimization of the Energy . . . . . . . . . . . . . . . . . . . . . 20
3.3 Variational Approach for Fermi Systems . . . . . . . . . . . . . . 22
3.4 HNC/0 Approximation For Fermions . . . . . . . . . . . . . . . . 23
vi
CONTENTS vii
4 Two and Three Dimensional Many Body Systems 26
4.1 3D Charged Boson Gas . . . . . . . . . . . . . . . . . . . . . . . . 26
4.1.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.1.2 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 27
4.1.3 Results and Analysis . . . . . . . . . . . . . . . . . . . . . 30
4.2 2D Charged Boson Fluid . . . . . . . . . . . . . . . . . . . . . . . 32
4.2.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.2.2 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 33
4.3 3D Spin Unpolarized Electron Gas . . . . . . . . . . . . . . . . . 36
4.3.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.3.2 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 37
4.4 2D Spin Unpolarized Electron Liquid . . . . . . . . . . . . . . . . 42
4.4.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.4.2 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 43
A Minimization of the Energy 50

CONTENTS viii
B The Equivalence of the Methods 54
C The Code 57
C.1 3D Electron Gas . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
C.2 3D Charged Bose Gas . . . . . . . . . . . . . . . . . . . . . . . . 59
C.3 2D Charged Bose Gas . . . . . . . . . . . . . . . . . . . . . . . . 61
C.4 2D Electron Gas . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

List of Figures
2.1 Radial distribution function for parallel and antiparallel spin con-
ditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.1 The representation of cluster functions . . . . . . . . . . . . . . . 15
3.2 The notations of classified cluster diagrams . . . . . . . . . . . . . 17
3.3 The diagrammatic representation of the operation between nodal

functions and convolution products . . . . . . . . . . . . . . . . . 18
3.4 An elementary diagram sample with four point structure. . . . . . 20
4.1 The plot of the convergence of the radial distribution function . . 29
4.2 Radial distribution function versus r/rs a for 3D Boson gas . . . . 30
4.3 Static structure factor versus kr0 for 3D Bose gas . . . . . . . . . 31
4.4 Radial distribution function versus r/rs a for 2D Boson gas in high
density regime. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.5 Radial distribution function versus r/rs a for 2D Boson gas . . . . 35
4.6 The plot of the convergence of the radial distribution function . . 39
ix
LIST OF FIGURES x
4.7 Radial distribution function of 3D electron gas . . . . . . . . . . . 41
4.8 The plot of the convergence of the radial distribution function for
2D electron gas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.9 Radial distribution function of 2D electron gas at several densities

for unpolarized 2D electron gas. . . . . . . . . . . . . . . . . . . . 46
4.10 Radial distribution function of 2D electron gas compared with

FHNC results for unpolarized case. . . . . . . . . . . . . . . . . . 46
List of Tables
4.1 Ground-state energy of the 3D charged boson fluid in Ry . . . . . 32
4.2 Ground-state energy results for 2D charged boson gas in Ry . . . 36
4.3 Correlation energy results for the unpolarized electron gas in Ry . 41
4.4 Correlation energy results for the unpolarized electron gas in high
density regime in Ry . . . . . . . . . . . . . . . . . . . . . . . . . 42
xi
Chapter 1
Introduction
After science has stepped to quantum mechanics from classical mechanics,

nanoscopic examination of matter has attracted extreme interest. The physi-
cal problems which include large number of interacting particles that is related
to quantum mechanics are named in general as quantum many body problems
[1]. The treated problems in this thesis are the 2D and 3D electron liquid sys-
tems. We have studied the radial distribution function and ground state energy
calculation of two and three dimensional electron liquids by the method of Kallio.
We additionally tried the method for two and three-dimensional boson systems
too.
The outline of the thesis is as follows. In the second chapter we concisely intro-
duce some background concepts to gain a clear perspective to the thesis: Jellium
Model, Wigner-Seits radius, particle density and Coulomb interaction are de-
fined, radial distribution function and static structure factor are discussed, and
an outline of the description of a quantum many body system is given to show
the systematic approach of the hypernetted chain method (HNC) that we have
used in this study.
In Chapter 3 we present the HNC Method aimed to investigate the distribution of

electrons and the energies of some simple models. We introduce how the method
1
is developed via the cluster diagrams in Sec. 3.1. Then, the founding of HNC/0
method, and its iteration process is exhibited in Sec. 3.2. Since, perturbation
theory only works for high electron densities (free electron gas), new methods
like coupled cluster method and variational methods were developed to solve the
problem of strong correlation [2]. One of them is the Fermi hypernetted chain
theory (FHNC) which is described in Sec. 3.3. In this method, one deals with
four highly nonlinear differential equations to find the ground state properties of
a Fermi liquid system. Also, it offers a hard and not exact way to evaluate the
elementary diagrams which are represented in Sec. 3.1. Therefore, in 1996, Kallio
and Piilo [3] offered a new technique to count the elementary diagrams in a much
simpler and exact way. This method, explained in Sec. 3.4, is basically using the
HNC/0 method which is modified for fermions to investigate the ground-state
properties of many-body fermionic systems.
In Chapter 4, the iterative schemes, the codes with explanations and the results of
the calculations for 2D and 3D electron systems are presented. We tried initially
for 3D charged bosons in Sec 4.1. Then, we moved our calculation to 2D charged
bosons in Sec 4.2. We defined the problems in detail, presented our results, and
then compared them with other results in the subsections of each section. In Sec
4.3, we produced the results for 3D free electron gas, and in the following section
we moved to 2D electron gas system to calculate the same properties with a fixed
scheme for their specific case.
2
Chapter 2
Background
The theoretical studies on electron liquids had an important role in the improve-
ments on novel materials and modern electronic devices [1]. Considering solid
state applications, it is preferred to solve the many body problems by assum-
ing the distribution of electrons is homogeneous. Such a model is the Jellium
Model defined in Sec. 2.1. In a many body system, the investigation of the
properties of the materials is easier by using the particle density instead of cal-
culating electron-electron interactions for each electron pair. For this reason, we
assume a homogenous electron gas with positive background, and investigate its
distribution and energy.
Our approximations are:

1) Three body correlations between electrons are neglected.
2) Jellium model is used.
3) Variational approach is used.
4) The Euler-Lagrange equation that we handle is a Schrödinger-like equation.
Hence, we assume that we have a Schrödinger equation with zero energy, and its
probability is the radial distribution function of the system.
5) Before the iterative calculation process, we presume that the Coulomb poten-
tial of the system is zero. This is to build up the fermionic correction term in the
Schrödinger-like equation.
3
First, I would like to introduce the crucial concepts that are needed to be under-
stood before dealing with the interacting electron gas problems.
2.1 Jellium Model
The Jellium Model is known as the homogenous electron gas with a uniform
positive background. Scientists are able to study the interaction effects by this
model, since it eliminates the effects coming from the periodic crystal potential.
In order to reach this homogenous case experimentally, scientists prepare a clear
and smooth environment in which electrons move freely and interact only among
themselves. Hence, the analytic and experimental results can be compared in this
model.
2.1.1 Particle Density and Wigner-Seitz Radius
Many elemental metals and semiconductors can be used for 2D and 3D electron
liquid studies. Metals are an example of 3D case since their electrons are able to
move freely, due to the almost spherical surface of conducting electrons which is
called the Fermi sphere. Every metal has a specific electron density which can
be defined by the number rs . The Wigner-Seitz radius (rs ) is the dimensionless
radius of the sphere that encloses one electron if we take average in a homogenous
environment. The rs value multiplied by Bohr radius gives the actual radius of
the imaginary sphere of one electron in units of Å. Bohr radius can be defined as:
~2
a0 = = 0.529Å. (2.1)
me2
The metallic rs range is roughly between 1.8 to 6 [1, p. 7-8]. However, by us-
ing doped metals one can change the electron density of the metal. To control
the electron density, we can use dopant concentration as a variable. In those
e2 e2
systems the Coulomb interaction formula changes from r
into r
where is the
4
static dielectric constant. Experimentalists prefer pure semiconductors due to
the disorder problems caused by dopant atoms. Pure semiconductors are at zero
temperature.
Density of particles in the system can be expressed in terms of rs and Bohr radius
(a0 ) for one, two and three dimensional systems as follows [1, p. 17] :

4π


 3
(rs a0 )3 3D
1 
= π(rs a0 )2 2D (2.2)
n 

2r a
s 0 1D.
2.2 Coulomb Interaction
e2
The Coulomb interaction between two electrons is vc (r) = r
in Gaussian units
for both 2D and 3D systems. To define the Coulomb potential in k space, one
takes its Fourier transform, and the result is as follows:

 4πe2 2 3D
q
vq = (2.3)
 2πe2 2D.
q
In the very high density range such that rs → 0, Coulomb potential between
particles can be neglected due to the dominance of the kinetic energy. We use
this fact as a first approximation in our calculations of highly dense systems. By
this way, we are able to set a clear analytic method, explained in Chapter 3,
to solve for the ground state energy and radial distribution function. Then, we
add the contribution of the Coulomb potential to the induced interaction of the
system so that the results become reasonable in the rs range from 0 to 20. It is
a useful density range for most solid state applications.
5
2.3 Radial Distribution Function
Radial distribution function describes the relative behavior of particles, and can
be defined as the distribution of electrons about any electron if the average is
taken over the motion of the particles [4, p. 60]. It gives significant clues about
the system studied. For instance, in order to calculate the ground state energy of
a many particle system, we need to obtain its exact radial distribution function.
Radial distribution function takes into account the two, three, and more particle
correlations. For the high density systems, we can neglect three and more body
interactions. The radial distribution function for only two body correlations is
called pair distribution function. It is defined as the probability per unit volume
that an electron is at r given that there is one at r = 0. The spin of the electron
has crucial importance on the g(r) of the system. There are two types of g(r)
if the electron gas is spin-unpolarized: g (r) = g (r) and g↑↓ (r) = g↓↑ (r). For
instance, g↑↓ (r) means the probability per volume such that the electron at r
has spin-up whereas the center electron has spin-down. These pair distribution
functions give the average state of the electrons. In order to show the calculation
of the radial distribution function, we start with a Slater Determinant:

φλ1 (r1 ) φλ1 (r2 ) ... φλ1 (rN )

φλ (r1 ) φλ (r2 ) . . . φλ2 (rN )
2 2
ψλ1 ,λ2 ...λN (r1 , r2 . . . rN ) = . .. .. ..
..

. . .

φλN (r1 ) φλN (r2 ) . . . φλN (rN )

where λj stands for jth state. Its square gives the N -particle density matrix.
However, we need only two particle density matrix to define the pair distribution
function of N particles. Therefore, it is found by integrating N -particle density
over all in two space coordinates:
Z
2
gss0 (r1 , r2 ) = v d3 r3 d3 r4 . . . d3 rN |ψλ1 ,λ2 ...λN (r1 , r2 . . . rN )|2 , (2.4)
where v is the fixed volume [4, p. 301]. We can sum over all possible distribution
6
wave functions by assuming that one-electron orbitals(φλj ) are orthogonal:
2
v2 X φλi (r1 ) φλi (r2 )
gss0 (r1 , r2 ) = .
N (N − 1) λ ,λ φλj (r1 ) φλj (r2 )

i j
We can define the one-electron functions as plane waves:

χs
φλ (r) = √ eik·r , (2.5)
v
where χs stands for the spin functions. Since the radial distribution function
defines the orbital section, the average of the spin functions is taken as: hχ↑ χ↑ i =
1, hχ↑ χ↓ i = 0, etc. For this two spin conditions, g(r) can be calculated as:
1 X
g↑↓ (r1 − r2 ) = (|eik1 ·r1 +ik2 ·r2 |2 + |eik1 ·r2 +ik2 ·r1 |2 )
N (N − 1) k ,k
1 2
2 (2.6)
2 N 1
= = ,
N (N − 1) 2 2
1 X
g (r1 − r2 ) = (|eik1 ·r1 +ik2 ·r2 + eik1 ·r2 +ik2 ·r1 |2 )
N (N − 1) k ,k
1 2
2 X
= [1 − ei(k1 −k2 )·(r1 −r2 ) ] (2.7)
N (N − 1) k ,k
1 2
1
= [1 − Λ(r1 − r2 )2 ],
2
2 X
where Λ(r) = [1 − eik·r ]. (2.8)
N k
The sum in the Eq. (2.8) is calculated over occupied states. As it is shown in
Eq. (2.6), the pair distribution function for the antiparallel spin case (g↑↓ ≡ g↓↑ )
equals to 12 . In the parallel spin case, λ(r) can be calculated as:
Z ik·r Z kF
2 3 e 3 3
λ(r) = dk n k = 3 ksin(kr)dk = j1 (rkF ). (2.9)
n0 (2π)3 rkF 0 rkF
One can see the pair distribution functions for parallel and antiparallel spins in
the Fig. 2.1. Therefore, the total radial distribution function which varies from
0 to 1 can be expressed as the sum of these two conditions:
g(r) = g (r) + g↑↓ (r). (2.10)
7
Figure 2.1: The pair distribution functions for parallel and antiparallel spin con-
ditions are found by Hartree-Fock approximation. The blue line represents the
g(r) for parallel spins, and the red line represents g(r) for antiparallel spins. The
functions are defined in [4, p. 309].
The definition is the same if the center electron is spin-down.
The pair distribution function has two remarkable features. First one is the
normalization rule of g(r):
Z
n0 d3 r[g(r) − 1] = −1. (2.11)
Second is one can calculate the Coulomb energy if the g(r) of the system is known:
e 2 n0 g(r) − 1
Z
EC = d3 r . (2.12)
2 r
2.3.1 Static Structure Factor
Structure factor can be found experimentally by the examination of scattering of

the x-rays, and electron or neutron beams from the semiconductor crystals. Thus,
the radial distribution function can be determined numerically by the relation
8
between S(k) and g(r):
e−ik·r
Z
1
g(r) − 1 = d3 k [S(k) − 1],
n (2π)3
Z (2.13)
3 −ik·r
S(k) − 1 = n d re [g(r) − 1].
Likewise the radial distribution function, static structure factor (S(k)) approaches
unity as k becomes large.
2.4 Description of a Quantum Many Body Sys-

tem
It is challenging to solve the many body Schrödinger Equation in highly dense

systems. Due to having many particles in a small area or volume, there will be
too many matrix elements so one cannot use the regular perturbative methods
and solve them. For this reason, high correlations in short range interactions
must be inserted into the wave function initially. We can denote the trial wave
function as ψT which is used for N particles in the system:
ψT (1 . . . N ) = F (1 . . . N )φ(1 . . . N ), (2.14)
where F is the operator factor which represents correlations. It is symmetric

under the exchange of particles, and φ is the totally symmetric wave function for
the free gas that does not take into account the interactions. The ground state
means all particles are in their lowest energy state such that φ(1 . . . N ) = Ω−N/2 .
The function φ is considered as the Slater determinant for fermions:
φ(1 . . . N ) = det |ψαi (j)|. (2.15)
Since the function φ is for the noninteracting free gas, the single particle wave
function ψαi (j) is taken as plane wave:
1 ikα ·rj
ψαi (j) = e i σ(j), (2.16)
Ω1/2
where σ can be spin-up | ↑i or spin-down | ↓i. Hence, the degeneracy is 2 for
unpolarized systems, and 1 for polarized systems. In the thermodynamic limit,
9
the volume which determines the boundary conditions goes to infinity. Now that
the trial wave function is specified, it can be used to get more information about
the many particle solid state system. In the variational approach, initially, a
trial wave function must be chosen. Through calculating the expectation value
of the Hamiltonian of the system, the upper bound to the ground state energy
can be found. The following variational method (Hypernetted Chain/0 Method)
to evaluate the ground state energy holds for a bosonic many particle system.
Since the wave function of a system is symmetric for bosons and antisymmetric
for fermions, the calculations are much easier for the bosonic case. The trial
wave function must be chosen carefully in order to use the variational approach
efficiently. Therefore, the Jastrow correlation wave function (f2 ) is selected for
two body correlations for the following reasons: It is zero or nearly zero if the
distance between two particles is less than the interval of the repellent part of
the potential. It approaches unity when the distance is large enough, to take no
account of correlation effect. Then, the free parameters of the Jastrow correlation
function can be found by minimization of the energy so that the optimum trial
function can be established:
δ hΨJ |H|ΨJ i
= 0. (2.17)
δf hΨJ |ΨJ i
Here, one can find the expectation value of the Hamiltonian through Monte Carlo
methods by handling the integrals in the energy equation directly. Monte Carlo
is a numeric method which gives exact results for the trial wave function. For this
reason, we will use the recorded results of Monte Carlo for our trial systems as a
check. In our calculations we use hypernetted-chain approximation, a more ana-
lytical method. Before founding the trial wave function, we should calculate the
expectation value of the Hamiltonian of the many body system. For this, we per-
form the cluster expansion of the integrals and sum up the terms systematically
by the aid of diagrammatic notation. Initially, we can calculate the expectation
value of the potential energy:
P
hV i 1 hΨJ | i<j V (rij )|ΨJ i N (N − 1) hΨJ |V (r12 )|ΨJ i
= = . (2.18)
N N hΨJ |ΨJ i 2N hΨJ |ΨJ i
10
Since we have chosen two particles out of N particles, we have to perform the
integrals over r3 , r4 . . . rN :
N − 1 dr1 . . . drN Ψ∗J (r1 , r2 . . . rN )ΨJ (r1 , . . . rN )V (r12 )

R
hV i
= R
N 2 dR|Ψ|2
R (2.19)
n2 N (N − 1) dR12 |ΨJ |2
Z
= dr1 dr2 V (r12 ) R ,
2N n2 dR|ΨJ |2
where dR12 = dr3 dr4 . . . drN . In Eq. (2.19) we notice that the radial distribution
function coincides with the second term in the integral:
R
N (N − 1) dR12 |ΨJ (r1 , . . . , rN )|2
g(r1 − r2 ) = R . (2.20)
n2 dR|ΨJ (r1 , . . . , rN )|2
Hence, we can display the potential energy as:

hV i
Z
n
= drV (r)g(r). (2.21)
N 2
The expectation value of the kinetic energy is found similarly:
N
hT i 1 X ~2 hΨJ |∇i 2 |ΨJ i ~2 hΨJ |∇1 2 |ΨJ i
=− =− . (2.22)
N N i=1 2m hΨJ |ΨJ i 2m hΨJ |ΨJ i
It is calculated by A. Polls and F. Mazzanti through the Jackson Feenberg Iden-

tity. The detailed solution is in [2, p. 59,60], and the final result is
hT i n ~2
Z
=− drg(r)∇2r ln f (r), (2.23)
N 2 2m
where f (r) is the Jastrow function. Consequently, the expectation value of the
Hamiltonian per particle of a many body system can be written as:
~2 2
Z
E n
e= = drg(r) V (r) − ∇ lnf (r) . (2.24)
N 2 2m r
The ground state wavefunction of a Bose system can be defined in the exponential
form:
1
Ψ0 (r1 , . . . , rN ) = exp χ(r1 , . . . , rN ) . (2.25)
2
In the Eq. (2.25), χ can be decomposed into its components as:
N
X X
χR (r1 , . . . , rN ) = un (ri1 , . . . , rin ), (2.26)
n=1 i1 ,...,in
11
where n is the number of interacting bodies, and u is the function of them. By
denoting 2.26 into 2.25 the Feenberg ansatz equation can be found:
" N #
1X 1 X
Ψ0 (r1 , . . . , rN ) = exp u2 (rij ) + u3 (ri , rj , rk ) + . . . . (2.27)
2 i<j 2 i<j<k
The equation can also be represented as:

Y Y
Ψ0 (r1 , . . . , rN ) = f2 (rij ) f2 (ri , rj , rk . . . ), (2.28)
i<j i<j<k
where fn = eun /2 . Now, since we assume that we are in a highly dense system of
particles, the many body interactions which include three and more particles can
be ignored. Thus, it is reasonable to choose the following Jastrow wave function
which considers only two-body correlations as the first guess of the bosonic wave
function:
Y
ΨJ (r1 , . . . , rN ) = f2 (rij ). (2.29)
i<j
As stated before, the Jastrow function approaches unity as the distance between
the two bodies approaches infinity due to the cluster property.
f (2) (r → ∞) → 1. (2.30)
12
Chapter 3
Hypernetted Chain Method
3.1 Definition
The hypernetted-chain (HNC) method is a variational method which was de-

veloped as an alternative to the field theoretical approach to cope with strong
correlations in many-body systems. It offers an analytical and systematic way
to improve the ground-state wave function [3] from first principles, and gives
a satisfactory description for Bose and Fermi systems with strong interactions,
where the perturbation theoretical method often fails. The HNC method is first
developed in 1959 by Leeuwen et al. [5] for classical liquid systems to understand
the thermo fluctuations of many particle systems at finite temperature, and it is
studied for classical systems mostly in 60s and 70s. After the developments in
quantum mechanics, HNC theory was also applied to quantum liquid systems to
understand quantum fluctuations at zero temperature in highly correlated Boson
systems such as Helium. In this thesis, we will take only the quantum applications
of this theory into consideration. The method works accurately for homogenous
quantum many-body fluids in the thermodynamic limit at 0 K. By the aid of
the HNC method, one can find the ground-state properties of interacting sys-
tems, such as ground-state energy, the radial distribution function, and the static
structure factor. An interacting system includes correlations between two, three
13
or more particles, hence there occur multidimensional integrals in the definition of
radial distribution function. The multidimensional integrals are known as cluster
functions. HNC uses a special schematic representation technique to deal with
these type of cluster integrals and makes the problem analytically solvable.
In Sec. 2.4, we showed that the radial distribution function is necessary to cal-
culate the total energy of a many particle system as summarized by Eq. (2.24).
Now, our aim is to clarify the description of the radial distribution function which
is given in Eq. (2.20) by defining the many-particle wave function in terms of the
Jastrow function as in Eq (2.29). Thereby, the pair distribution function can be
written as follows:
RQ
N (N − 1) f 2 (rij )dR12
g(r) = R Qi<j 2 2 . (3.1)
n2 i<j f2 (rij )dR
In order to solve the equation easily, we will introduce the cluster function h(r) =
f22 − 1 which goes to zero as the distance between particles approaches infinity.
RQ
N (N − 1) (h(r) + 1)dR12
g(r) = R Qi<j . (3.2)
n2 i<j (h(r) + 1)dR
After applying power expansion, g(r) becomes:

R P P
N (N − 1)Ω−N dR12 (1 + i<j h(rij ) + i<j,k<l h(rij )h(rkl ) + · · ·)
g(r) = R ,
n2 Ω−N
P P
dR(1 + i<j h(rij ) + i<j,k<l h(rij )h(rkl ))
(3.3)
where Ω−N stands for the total volume of the system with N particles. The sum
can be converted to its integral form. Here is an example:
Z N
X Z
−N +2
dR12 Ω h(r1i )h(ri2 ) = n drj h(r1j )h(rj2 ). (3.4)
i=3
In Eq. (3.4), the reason that the sum starts from i = 3 to N is that we chose
two external particles (i = 1, 2) and look for the internal particle that integrates
these two external particles. Because the number of remaining particles is N −
2, the volume factor becomes Ω−(N −2) . We take the integral with respect to
dR12 = dr3 , dr4 , ..., drN , and sum up the multiplication of two h(r) functions for
N − 2 different points. This sum can be written as an integral simply by taking
infinitesimal points. At the right hand side of Eq. (3.4), n is the density of
14
particles, and the integral is taken with respect to rj , where j stands for the
internal particle.
The multidimensional integrals including one or more cluster functions are called
cluster terms (see Eq. 3.3). As the amount of the integrated coordinates in-
creases, the calculation of the cluster terms become more complicated. We will
introduce the diagrammatic notation of cluster integrals to clarify the problem.
This notation leads us to realize the integrals with the same result and the essen-
tial cancellations between the infinite sets of cluster functions so that the problem
becomes solvable. Below are the rules of the diagrammatic analogy:
• One diagram can represent many cluster integrals.
• Every internal (integrated) particle corresponds to a solid circle.
• Every external (not integrated) particle corresponds to an open circle.
• Every h(rij ) function connecting the particles i and j, is represented as

dashed line.
The difference in the integrated coordinates does not change the result of the
integral. For instance, in Eq. (3.4), the integral of h(r15 )h(r25 ) and h(r18 )h(r28 )
gives the same outcome, and they are delineated by the same diagram. We
can represent the integral in the Eq. (3.4) as in Fig. 3.1(a) where 1 and 2
are the particles which are integrated. Further, the internal point, j, does not
change the yield of the integral. The diagram notation of the cluster function
Figure 3.1: The representation of cluster functions. a) for two integrated particles,
b) for two nonintegrated particles, reproduced from [2, p. 64].
15
h(r12 ) (N −2)(N −3) R
2Ω
drij h(rij ) is depicted in the Fig. 3.1(b) as an example. Here,
rij is the distance between particles i and j. In this cluster function, h(r12 ) is not
involved in the integral, since the coordinates 1 and 2 are not integrated. The
factor (N − 2)(N − 3) is multiplied due to the selection of two points from the
number of possible choices. Also, the term is divided by Ω that comes from the
R R R
change of variables, such that dr1 dr2 → dRdr = Ω dr where r = r1 − r2
and R = (r1 + r2 ). Additionally, the term is multiplied by 1/2 due to the fact
that the result is the same if the label of the particles i and j are exchanged.
The other diagrams can be evaluated in this way. We note that the particles
should be close enough to consider their correlation effect, since f (r → ∞) → 1
corresponds to h(r → ∞) → 0. This means the diagrams including more lines
contribute less to the radial distribution function. Therefore, we can use the less
complicated diagrams to approximate g(r). The diagrams can be classified as
linked or unlinked. Linked diagrams have only one way to connect the particles
which are represented as open or closed circles in the diagram. However, there
are some not connected parts in the unlinked diagrams. The unlinked and linked
diagrams can be seen in the Fig. 3.2. Linked diagrams are grouped as reducible
and irreducible. In reducible diagrams there is an articulation point which is a
circle pointed with an arrow. It separates the diagram into different parts [2,
p. 65]. It is hard to decide which of those diagrams should be kept for the radial
distribution. Thankfully, there is a theorem which states that we can consider
only the cluster functions corresponding to irreducible diagrams in the numerator
of Eq. (3.3), since the reducible and unlinked diagrams on the numerator are
cancelled by the denominator up to the order 1/N [6]. Thus, one can simply define
the pair distribution function as the sum of topologically distinct irreducible
diagrams, as in Eq. (3.12). The irreducible diagrams are divided into composite
diagrams labeled as X, and simple diagrams which are separated to nodal and
elementary diagrams labeled as N and E, respectively. In a diagram, if there is
only one line from one external point to another, then the internal points on this
line are denoted as nodes. In Fig. 3.2, it can be seen that elementary diagrams
have no node whereas nodal diagrams have at least one node.
16
Figure 3.2: The notations of linked, unlinked, reducible, and irreducible diagrams,
reproduced from [2, p. 66].
The cluster diagrams are formed by the aid of convolution and algebraic prod-
ucts. The integrals of functions a(rik ) and b(rkj ) are represented as convolution
products: Z
n(a(rik )|b(rkj )) ≡ n drk a(rik )b(rkj ), (3.5)
where k is the nodal point, because it is the only connection between the two
external points: i and j. In addition to this, the algebraic product of the two
functions produces a composite i − j subdiagram. Herein, we can form the HNC
equation, but omit the contribution of elementary diagrams initially. N (r) is the
sum of the nodal diagrams, and we take solely the subset of them found by the first
iteration. In order to build the chain, first the convolution product of N (r) and
h(r) is carried out. Thus, we obtain all nodal diagrams with the exception of the
convolution product of two h(r) functions [2, p. 68]. This situation is exhibited in
the Fig. 3.3. Hence, the selected set of diagrams satisfies the following equation:
17
Figure 3.3: The diagrammatic representation of the operation between nodal
functions and convolution products, reproduced from [2, 66].
n(h(rij )|N (rj2 )) = N (r12 ) − n(h(rij )|h(rj2 )). (3.6)
By this way, N (r) is constructed, and it is used to establish composite diagrams.

Composite diagrams include at least two 1 − 2 subdiagrams. The sum of these
N n (r)
subdiagrams consists of n!
terms where n is the number of 1 − 2 subdiagrams.
The terms are divided by n! in order not to count the same terms due to the
symmetry. Hence; in the thermodynamic limit, the sum can be written as
∞
X Nn
= eN − N (r) − 1. (3.7)
n=2
n!
If we multiply exp[N (r)] with h(r) to obtain more composite diagrams, we get
the same diagrams except the term eN h(r). Thus, the sum of all the composite
diagrams which has the notation X(r) equals:
X(r) = eN − N − 1 + (eN )h(r). (3.8)
Now, we can place the equation h(r) = f 2 (r) − 1 into Eq. (3.8), and find X(r)
as:
X(r) = f2 2 (r)eN (r) − N (r) − 1. (3.9)
Here, we see the results of the first iteration in the HNC scheme: X(r) and N (r).
We place the newly found X(r) into the Eq. (3.8) instead of h(rij ):
n(X(r1j )|N (rj2 )) = N (r12 ) − n(X(r1j )|X(rj2 )). (3.10)
18
By this equation, N (r) can be found directly, and used in Eq. (3.9) to find
the new X(r). This is the basic HNC iteration process between equations (3.8)
and (3.9) that lasts until the convergence is reached. Since we did not take the
elementary diagrams, it is named as HNC/0 method. However, now we can add
the contribution of elementary diagrams to the term X(r):
X(r) = f22 (r)e[N (r)+E(r)] − N (r) − 1. (3.11)
As stated before, g(r) is the total sum of the topologically distinct diagrams which
titled as nodal, elementary and composite diagrams.
g(r) = 1 + X(r) + N (r) = f22 (r)eN (r)+E(r) . (3.12)
The iteration scheme for the HNC Method accounting elementary diagrams is
generated by the equations (3.10) and (3.11). Since Fourier transform converts
convolutions to algebraic products, it aids to solve N (r) in Eq. (3.10).
X 2 (k)
N (k) = . (3.13)
1 − X(k)
At this point, we will define the Fourier transform as:

Z
a(k) = n dreik·r a(r), (3.14)
where n is the particle density. Static structure factor, defined in Sec. 2.3.1, can
also be expressed in terms of cluster functions by using Eq. (3.11) due to the
Fourier relationship described in Eq. (2.13):
S(k) − 1 = F[g(r) − 1]
= F[(1 + X(r) + N (r)) − 1] (3.15)
= X(k) + N (k).
Now, we can define the functions N and X in the k space in terms of the static
structure factor:
(S(k) − 1)2
N (k) = ,
S(k)
(3.16)
S(k) − 1
X(k) = ,
S(k)
19
Figure 3.4: An elementary diagram sample with four point structure.
where X(r) may not include the contribution of elementary diagrams depending
on the HNC process. The elementary diagrams consist of at least four points,
and the simplest one is represented in the Fig 3.4. It is not possible to express
sum of the infinite set of elementary diagrams in a closed form. At this point,
one should consider making an approximation. One approximation is counting no
elementary diagrams and this approach is named HNC/0, as mentioned before.
On the other hand, one can consider only the four points elementary diagram
(see Fig 3.4) denoted as E4 , and this approach is labeled as HNC/4. However; if
we use g(r) − 1 in place of h(r) in the HNC/4 approach, the elementary diagram
contributions can be built on E4 after each iteration, because g(r) consists of the
sum of all the diagrams. Hereby, we can do a better approximation to the exact
result. Also, we note that HNC/0 method, which is iterated on g(r) − 1 instead
of h(r), gives a good approximation for the bosonic case in the high density limit.
Nevertheless; in the fermion case, the contribution coming from the elementary
diagrams is not negligible. This case will be treated in the Sec. 3.3.
3.2 Minimization of the Energy
As defined before in Eq. (2.24), the energy per bosonic particle depends on g(r)
and f2 (r) :
~2 2
Z
n
e= drg(r) V (r) − ∇ lnf (r) . (3.17)
2 2m r
Here, V (r) equals the Coulomb potential for the systems we studied. To solve this
equation, the radial distribution function and Jastrow function can be defined in
terms of each other if we remember Eq. (3.11).
g(r) = 1 + X(r) + N (r) = f22 (r)eN (r)+E(r) . (3.18)
20
Here, we can take the functional derivative of the energy with respect to g(r) or
f (r) to minimize the energy and find the optimum trial function. We thought
that it is better to get rid of f (r) term by defining it in terms of g(r) so that
we do not need to find the Jastrow correlation function. Instead, we will use
a reasonable g(r) as a first guess for the studied many body system. In this
manner, we can define our energy per particle in terms of g(r). The detailed
solution of the following process is expressed in Appendix A. First, we do the
HNC/0 approximation to solve the energy minimization problem handily, that
means E(r) is taken as zero in Eq. (3.11). After Fourier transforming N (k),
which is defined in Eq. (3.13), and taking the log of both sides we get:
(S(k) − 1)2 ik·r
Z
1 1
ln f (r) = ln g(r) − e dk . (3.19)
2 (2π)3 n S(k)
Then, we insert Eq. (3.19) into Eq. (3.17), and find the total energy per particle
as:
~2 − 1)3 n ~2
Z Z Z
n 1 2 (S(k)
e= drg(r)vc (r) − dkk − drg(r)∇2 ln g(r).
2 8m (2π)3 n S(k) 2 4m
(3.20)
p
It is easier to do the variational derivative with respect to g(r) instead of
taking the functional derivative with respect to g(r). Here, I should note that the
p
condition for optimization is equivalent for both trials. g(r) can be denoted as
G so that the energy minimization with respect to G is:
δe(n)
0=
δG Z
~2 3
~2

2 (S(k) − 1)
Z Z
δ n 2 1 2
= drG (r)vc (r) − dkk + n dr(∇G(r)) .
δG 2 8m (2π)3 n S(k) 2m
(3.21)
In the Appendix A, we have taken the functional derivative step by step by using
the Parseval identity. Then, we find the so called zero energy Schrödinger like
equation with the wave function G(r):
δ ~2 2
(e(n)) = 0 ⇒ − ∇ G(r) + [vc (r) + WB (r)]G(r) = 0. (3.22)
δG m
Here, WB (r) is the bosonic induced potential which is defined in k space as:
2
~2 k 2 (S(k) − 1)

WB (k) = − (2S(k) + 1). (3.23)
4m S(k)
21
In the energy equation we have the static structure factor S(k). Since we can
study the energy minimization in the momentum space, we can do the functional
derivative with respect to S(k):
δe(n) 1
=0 ⇒ S(k) = q , (3.24)
δS(k) 2
V (k)
(k) ph
where (k) is the kinetic energy, and Vph (k) is the Fourier transform of the particle-
hole effective interaction which can be represented as:
~2 p
Vph (r) = g(r)vc (r) + WB (r)[g(r) − 1] + |∇ g(r)|2 . (3.25)
m
The iteration process to find the optimum g(r) is summed up by the above equa-
tions (3.22 - 3.25). At first, one should use a reasonable first guess for g(r), and
corresponding S(k). First, we place S(k) into the Eq. (3.23) and Fourier trans-
form it to find WB (r). Then, we put g(r) and WB (r) into the Eq. (3.25), so that
Vph (r) is calculated. By Fourier transformation, one can find Vph (k) which will
be put into the Eq (3.24) so that the new S(k) is found. The iteration process
should continue until there will be no observable change in g(r). Then, the energy
per particle can be obtained by replacing the last g(r) and S(k) into the energy
equation for a bosonic system.
3.3 Variational Approach for Fermi Systems
If the many particle system consists of fermions, the contribution coming from
the elementary diagrams become so significant that we cannot eliminate them as
we do in the HNC/0 approximation. In the Fermi Hypernetted Chain Method
(FHNC), one must take at least four point elementary diagrams. After several
steps, we will introduce only the compact definition of the radial distribution
22
function which is obtained in detail in [2, p 102-111].
g(r) =f22 (r) exp[Ndd (r) + Edd (r)]

l2 (rkF )
[1 − + Nee (r) + Eee (r) + [Nde (r) + Ede (r)]2
v (3.26)
+ 2[Nde (r) + Eee (r)]v[Ncc (r) + Ecc (r)]2
+ 2l(rkF )[Ncc (r) + Ecc (r)]].
Here, l(x) is the Slater function, v is the spin degeneracy, and a, b, c, d notations
stand for the four points of E4 . A sample diagram of four point elementary
function was represented in Fig. 3.4 in which the dashed lines corresponds to
the cluster functions in the integral. Because this integral is too complicated to
be evaluated exactly, and easily; a more effective method is required. Hence, in
1996, Kallio and Piilo [3] have introduced a novel technique which offers an exact
way to count elementary diagrams for fermions. In this technique, we can use the
simple HNC/0 method with an additional correction term that accounts for the
contribution of elementary diagrams. As well as it offers an exact way to count
elementary functions, it reduces the time of the computation.
3.4 HNC/0 Approximation For Fermions
The HNC/0 method for fermions has developed in 1995 by Kallio et. al [7].
The method basically assume that there is a zero energy Schrödinger equation
p
which is satisfied by the probability amplitude ψ = g(r). The potential part
of the Schrödinger equation consists of the Coulomb potential (vc (r)) and the
induced interaction potential denoted as W (r). The induced interaction, which
is explained in Eq. (3.22) for the bosonic case, includes one additional term,
WF (r), for the fermionic case. It is the key part of the novelty of this method,
since WF (r) is the correction term to account the elementary diagrams exactly.
This new term is found analytically by using the fact that Coulomb potential
can be neglected at an extreme electron density situation. By the aid of the
fermionic correction term, we are able to use the HNC/0 method for fermions
such that the hardness and the time spent are thoroughly low if compared to the
23
FHNC method. At first, we will introduce the iteration process of Kallio that is
solved for the bosonic case, and then move the discussion to fermions. In this
way, we will use a defect function R to solve the highly nonlinear zero energy
equation (3.22). The iteration process is described in detail as follows: First, we
will define the most realistic radial distribution function for the handled many
body problem. Here, note that the wave function of the zero energy Schrödinger
p
Equation is equal to g(r). Then, the pair distribution function can be put into
the following defect function:
∇2 (g − 1) ∇2 ψ
R(r) = − + . (3.27)
2 ψ
By this way, the new S(k) can be found via placing the Fourier transformed R(r)
into the following equation:
1
S(k) = q (3.28)
4γ 4
1+ k4
− k2
R(k)
with γ is 4πn/a0 . a0 is the Bohr radius defined in Eq. (2.1) . Next, the new g(r)
is obtained through the relation between S(k) and g(r) if we recall Eq. (2.13):
e−ik·r
Z
1
g(r) − 1 = d3 k [S(k) − 1] = F −1 [S(k) − 1]. (3.29)
n (2π)3
This is the iteration scheme for bosons that is introduced in Sec. 3.2. Actually,
this is the same HNC/0 process, but there is a little difference in their ways to
solve the nonlinear zero energy Schrödinger equation and we show their equiv-
alence in Appendix B. Now, we can move our discussion to many-body fermion
case. At this point, we add a fermionic correction term to the induced potential
WB (r) in the Eq. (3.22):
~2 2
− ∇ ψ + (vc + WB + WF )ψ = 0. (3.30)
m
As mentioned before, WF (r) can be found by the assumption that we deal with an
extremely dense fermion system so that the Coulomb potential can be neglected:
~2 2
− ∇ ψ + (WB + WF )ψ = 0. (3.31)
m
Hence, the WF is able to be found directly by putting WB into the Eq. (3.31):
2
~2 k 2 (SF (k) − 1)
2 2
~ ∇ψ
WF (k) = − (2SF (k) + 1) + F , (3.32)
4m SF (k) m ψ
24
√
where ψ = gF . This time new S(k) is calculated in the form
1
S(k) = p . (3.33)
1 + 4γ/k + (4/k ) [(m/~2 )WF (k) − R(k)]
4 2
The iteration process is applied between Eq. (3.27) and Eq. (3.33). After the
convergence is reached, one may consider calculating the energy of the system
per fermion. The potential energy of a many particle system is calculated in
two ways. One may use g(r) or S(k) and get the same result for the interaction
energy:
Z
V n
= d3 r[g(r) − 1]vc (r)
N 2
(3.34)
4πe2
Z
1 3
≡ d k[S(k) − 1] .
2(2π)3 k2
In our calculations, we used both methods to compare and verify our energy
values.
25
Chapter 4
Two and Three Dimensional

Many Body Systems
We have investigated the ground state properties of charged bosons and electrons
including the calculation of ground state energy, radial distribution function, and
static structure factor. Bosons are particles with integer spin values, and they
obey Bose-Einstein statistics. As a bosonic system we used composite boson sys-
tem properties, and investigated the ground state behavior of 3D and 2D charged
bosons. Fermions obey Fermi-Dirac statistics and have half-integer spin, and they
are restricted by Pauli exclusion principle. Fermions include electrons, protons
and neutrons. As fermions, we studied 3D, and 2D electron liquid systems.
4.1 3D Charged Boson Gas
4.1.1 Definition
The charged boson fluid has been usually treated as a basic model of quantum
many-body systems in statistical mechanics. Even though the Bose fluid has no
26
applications in the real world, it has attracted much interest as a model for su-
perconductors, after the discovery of high temperature superconductivity. There
is also a possibility to experimentally generate them using metals such as palla-
dium or vanadium highly doped with deuterium [8]. Additionally, this model has
a relevance to astrophysics [9, 10, 11] in the expression of pressure-ionized helium
in cold degenerate stellar matter.
We first tried our method on many body charged boson fluid, since it is the
fundamental quantum system to which the hypernetted-chain method can be
applied. Charged boson gas is a many particle system which has only the Coulomb
interaction between its bosonic particles. We applied the main hypernetted-
chain method to calculate the pair distribution function, static structure factor,
and ground state energy of the bosonic system. Also, we investigated how the
distribution function and energy depends on rs . Additionally, we compared our
results with Diffusion Monte Carlo (DMC) and HNC results.
4.1.2 The Algorithm
For the three dimensional Bose gas, it is known that the radial distribution func-
tion equals to 1 for the noninteracting case. Since we do not include the Coulomb
potential which is the only interaction between charged boson particles; the prob-
ability to find another boson at any distance from our reference particle equals to
1. Therefore, gB (r) = 1. Due to the Fourier transform relationship between g(r)
and S(k) in Eq. (3.29), the corresponding static structure factor of that radial
distribution function will be 1 too. Here is the step by step description of the
iterative scheme of 3D Bose gas which is summarized in Sec. 3.2, and explained
in Sec. 3.4:
1 Initially, we use dimensionless variables to make our calculations simpler:
r/r0 → r, k · r0 → k (4.1)
where r0 = rs · a is the radius of a sphere that encloses one particle on
27
average. The defect function is defined as:
∇2 (g − 1) ∇2 ψ
R(r) = − + , (4.2)
2 ψ
p
where ψ = g(r). To find the defect function, we take the Laplacian of
g(r) in the Eq. (4.2):
1 ∂ 2 ∂g
∆g = (r ) (4.3)
r2 ∂r ∂r
We write rs = 0 in the code for the first calculation, since we are starting
with the noninteracting case or high density limit.
The defect function can be calculated by using gB (r) = 1, and SB (k) = 1.
After one iteration, we will use the new g(r) and S(k) to calculate R(r).
2 One can find the new S(k) through taking the Fourier transform of R(r),
and placing it into the equation:
1
S(k) = q (4.4)
4γ 4
1+ k4
− k2
R(k)
where γ = 4πn/a. Here, a is the Bohr radius defined in Eq. (2.1), and n is
the density of bosons:
3 1
n= (4.5)
4π (rs a)3
3 By inserting new S(k) function into the following equation, one can find
new g(r):
e−ik·r
Z
−1 1
g(r) − 1 = F [S(k) − 1] = d3 k [S(k) − 1]
n (2π)3
k 2 (S − 1) −1
Z Z Z 2π
4π 3 ikr cos θ
= − (rs a) dk d(cos θ)e dφ
3 (2π)3 1 0 (4.6)
2(rs a)3
Z
sin(kr)
= dk[S(k) − 1]k 2
3π kr
Z
2 sin(kr)
≡ dk[S(k) − 1]k 2
3π kr
4 Now, we return to step 1, and insert the new definitions of g(r) and S(k) into
the equations (4.2), and (4.3). This whole procedure should be repeated
until the convergence is reached. Fig. 4.1 shows that the radial distribution
28
function converges to a certain form after 60 iterations for a 0.1 increase in
rs . Although normally, 10 to 20 iterations are enough to reach convergence,
the iteration number is 60 for this specific case. This is because we used a
mix of the new g(r) and the previous g(r). For the 3D boson case when rs
is between 0 to 1, the new radial distribution function can be taken as:
4 · gprev (r) + gnew (r)
g(r) = (4.7)
5
Mixing a function is a known numerical tool used in iterative algoritms
to reduce errors, and it can be adjusted for different rs values. We have
plotted the convergence scheme for all of the calculations to determine the
mixing rate by checking if the convergence is smooth. In Fig. 4.1, a dot
R
corresponds to the value of the integral: dr|gaf ter (r) − gbef ore (r)|. We
stopped the computation when this value become smaller than 1 × 10−7 .
We plotted the last g(r) and its corresponding S(k) functions for different
rs values in Fig. 4.2.
Figure 4.1: The plot of the convergence of the radial distribution function
5 After finding the interacting state g(r), we calculate the potential energy
of the system first:
Z
V n
= V (rs ) = d3 r(g(r) − 1)vc (r) (4.8)
N 2
where vc (r) = e2 /r. Another way of finding the potential energy is to apply
the integration in k space by using the interacting S(k):
4πe2
Z
1 3
V (rs ) = d k[S(k) − 1] (4.9)
2(2π)3 k2
29
From Eq. (4.9), it is obvious that potential energy depends on the value rs .
We used both ways (4.8, 4.9) to calculate V (r) to check our results. Now,
we can find the ground state energy by coupling constant integration, using
Eq. (4.8) : Z rs
1
EGS = 2 drs0 rs0 V (rs0 ) (4.10)
rs 0
In our code, the upper limit of the integration is rs as seen in Eq. (4.10).
For rs = 5, the potential energy is computed from rs = 0 to rs = 5 by 0.1
steps. Then, V (rs ) is interpolated in order to have a continuous function
that will be integrated with respect to the variable rs . Hence, we could
calculate the total energy for rs = 5. The same calculation is done for
different rs values between 0 − 20 by placing the aimed rs value into the
Eq. (4.10).
4.1.3 Results and Analysis
The plot of the radial distribution function as a function of r/rs a can be seen in
the Fig. 4.2.
Figure 4.2: Radial distribution function versus r/rs a for 3D Boson gas
Due to the definition of radial distribution function, and the fact that we take g(r)
30
as the probability amplitude of the Schrödinger-like equation; we can assume that
it acts like the probability of finding a particle at r distance from the reference
particle. As expected, Fig. 4.2 shows that g(r = 0) approaches to zero as rs
increases, since the probability of finding two particles at the same place decreases
with the decreasing density.
The corresponding static structure factor versus kr0 is plotted for different rs
values which is shown in the Fig. 4.3.
Figure 4.3: Static structure factor versus kr0 for 3D Bose gas
The static structure factor takes values from 0 to 1 for all of our calculations.
Its reason can be easily verified via Eq. (4.4). As exhibited in Eq. (4.9), this
function is also used to find the total energy. Since there is no Fermi level for
bosons, we do not include the Fermi energy part in our energy definition. After the
computations, we recorded the ground state energy results for different rs values,
and then compared with the results with other methods. In Table 4.1, one can
see that our results are very close to the exact DMC results, although we did not
take three body correlations into account. We also compared our results with the
results of Apaja and Halinen who uses the same HNC/0 Method including triplet
correlations. The energy comparison from Table 4.1 shows that there is not a
significant difference between our results and theirs even though our method is
much simpler. This part of the study is important, since it demonstrates that
our iteration method and code works very well, and we are ready to convert our
31
rs Egpres EgHN C EgDM C
0.1 -4.46523 -4.48857
1 -0.76618 -0.77675 -0.77664
2 -0.44841 -0.45203 -0.45192
3 -0.32523
5 -0.21432 -0.21657 -0.216420
10 -0.12135 -0.12144 -0.121353
20 -0.06659 -0.06664 -0.066639
Table 4.1: Ground state energies of the 3D charged boson fluid in Ry at several
densities. Egpres present calculation, DMC results from Ref. [12], HNC results
from Ref. [8].
scheme into other bosonic and fermionic cases. Moreover, we calculated a highly
nonlinear differential equation by an analytically based method via iterations,
and mathematical operations. The integrals and derivatives during the iteration
process are evaluated numerically.
4.2 2D Charged Boson Fluid
4.2.1 Definition
As mentioned Sec. 4.1, the fluid of 2D charged bosons has not been experimentally
realized. However, its basic model has importance for real world applications due
to being a 2D many particle system which has Coulomb potential between its
particles. Moreover, the recent studies on quantum Hall effect show that the
electron liquid in strong magnetic field acts like a quantum many-body system.
To understand the basic principles of the quantum Hall effect, one can initially
consider the ground state of a 2D charged boson in the absence of an external
magnetic field.
32
4.2.2 The Algorithm
In this part of the study, we have calculated the ground state energy of 2D
charged Bose gas. We have converted the iterative scheme of three-dimensional
system into two-dimensional including integrations, Fourier transforms, and the
definition of particle density. A short description of the iteration process can be
seen in the following:
1 Here, we use the same unitless system which holds for 3D Bose gas. For
the first iteration, gB (r) and SB (k) are equal to 1.
The particle density for the 2D Bose gas:
1
n= (4.11)
π(rs a)2
By using the equivalence of two HNC methods which is mentioned in Sec.

3.4 and solved in Appendix C; the two-dimensional γ is found as:
2πn 2
γ= = 2 3 (4.12)
a rs a
The defect function R(r) is defined in Eq. (4.2). While we are calculating
R(k), we have changed the definition of Laplacian and Fourier transform
from 3D to 2D in the code. The laplacian is taken in polar coordinates:
1 ∂ ∂g
∆g = (r ) (4.13)
r ∂r ∂r
2 Now, we can find the new S(k) by placing the functions γ and R(k) into
the equation:
1
S(k) = q (4.14)
4γ2d 4
1+ k3
− k2
R(k)
3 To find new g(r), one can use the Fourier transform relationship between
33
S(k) and g(r):
e−ik·r
Z
−1 1
g(r) − 1 = F [S(k) − 1] = d2 k [S(k) − 1]
n (2π)2
Z π
(S − 1)
Z
2
= π(rs a) dk k dθeikr cos θ
(2π)2 0 (4.15)
(rs a)2
Z
= dk[S(k) − 1]kJ0 [kr]
2
Z
1
≡ dk[S(k) − 1]kJ0 [kr]
2
4 As stated in Sec. 4.1.2, we use the mixed g(r) in step 1, and continue the
R
iteration process until the value of the integral: dr|gaf ter (r) − gbef ore (r)|
become less than 1 × 10−7 .
5 The potential energy can be found as follows:
e2
Z Z
n 2 1
V (rs ) = d r(g − 1)vc (r) = rdrdθ(g − 1)
2 2π(rs a)2 r
2 2
Z Z
e e 2
= 2
dr(g − 1) = dr(g − 1) (4.16)
(rs a) 2a rs 2 a
Z
2
≡ dr(g − 1) Ry
rs
It is defined in terms of Ry, the Rydberg constant, which equals to e2 /2a.

Then, we can find the total energy by the same way decribed in the step 5
in Sec. 4.1.2.
In Fig. 4.4 the radial distribution function versus r/rs a is plotted for different
boson densities. As stated previously, we initiated the iteration from the high
density limit rs → 0, where g(r) = 1, to the lower densities. When rs increases,
the contribution of the Coulomb interaction on the ground state energy also
increases. At the high density limit, the kinetic part of the energy expression
dominates over the potential part. Consequently, the Coulomb interaction be-
tween particles does not contribute a significant amount to the total energy of
34
Figure 4.4: Radial distribution function versus r/rs a for 2D Boson gas in high
density regime.
the system. That is why we can assume that Coulomb force is switched off for
very low rs values. Oppositely, when rs becomes higher, the potential part of
the energy expression dominates. Hence, the Coulomb potential become more
important at higher rs values as you can see in Fig. 4.5. Despite of the fact that
Figure 4.5: Radial distribution function versus r/rs a for 2D Boson gas
rs must be much higher than 20 in the case of a Wigner crystal, we can still see
the effect of the crystal state in Fig. 4.5. As we know, the plot of a Wigner
crystal must have many peaks in its g(r) plot. Nevertheless, since we are still
35
far away from that case, we only see a small peak on the function while we are
approaching to the crystal state. That is why we observe a higher peak for a
higher rs value in Fig. 4.5.
The numerical results for 2D bosons are given in Table 4.2. Actually, we can
rs Ecpres EcST W EcST LS

0.1 -5.8322
1 -1.1337 -1.1062 -0.1103
2 -0.6640 -0.6631 -0.6484
3 -0.4872 -0.4818 -0.5965
4 -0.3841 -0.3796
5 -0.3174 -0.3133 -0.3078
10 -0.1732 -0.16685 -0.1724
20 -0.0916 -0.086024
Table 4.2: Ground-state energy results for two-dimensional charged boson gas
in Ry. STW results are from Ref. [13] using a parametrized variational wave
function, STLS results from Ref. [14] calculated using the STLS scheme.
safely say that three-body correlations are more important for 2D systems than
3D systems due to the fact that the probability to find three particles close to
each other is higher for the particles which are restricted in two dimensions. Al-
though our results for the ground-state energies are not ”exact” due to neglecting
triplet correlations, they are very close to the other numerical data in the various
calculations which are summarized in Table 4.2.
4.3 3D Spin Unpolarized Electron Gas
4.3.1 Definition
Electron gas is a model consisting of electrons interacting by Coulomb forces.

Owing to the article published in 1957 by Gell-Mann and Brueckner [15] on the
correlation energy of electron gas at high densities, there has been a growing
36
interest to understand the properties of electron gas on positive uniform back-
ground. The studies on homogenous electron gas are particularly significant in
condensed matter physics [16, 17] due to the fact that electron gas can be used
as a simple reference model for electronic structure calculations.
As it is mentioned in Chapter 1, there are several methods to solve many electron

systems, and we handled FHNC which provides an analytical solution of a varia-
tional theory based on Jastrow function. In Chapter 3, we explained how FHNC
is developed, and the problems related to elementary diagrams. In Sec. 3.4, we
introduced the improvement on FHNC by Kallio et al. who propose a simpler
way to achieve more accurate results. In this work, our main aim is to contribute
to the understanding of the electron gas, and to generalize the usage of this simple
method to solve other Fermi systems such as two-dimensional electron liquids.
For this reason, we first reproduced the results found by Kallio, and then applied
this method for 2D electron gas. In Sec. 4.1.3, we had mentioned that DMC
method provides exact results for bosons, so that we could check our analytical
results against them. However; for fermions, DMC results are not exact because
of the Slater determinant nodes. Nevertheless, one can still compare their results
with other methods to evaluate the accuracy of their method.
4.3.2 The Algorithm
In Sec. 3.4, we have described the novel method of Kallio and Piilo, and explained
the trick they used to find the elementary diagrams between equations (3.30,
3.32). Here, we will show how this method is applied for unpolarized 3D electron
gas by discussing the iteration scheme step by step:
1 First, we use dimensionless variables to simplify our problem.
r/rs a → r, k/kF → k (4.17)

1
where kF is the Fermi wave vector. For 3D electron gas kF = ( 9π
4
) 3 rs1a . As
mentioned before, we write rs = 0 in the code for the starting case, and
37
write the known definitions of the noninteracting g(r) and S(k). For the
unpolarized electron gas; these functions are:
2
9 j1 (kF r)
gF (r) = 1 − . (4.18)
2 kF r

3 k 1 k 3
− if kkF ≤ 2

4 kF 16 kF
SF (k) = (4.19)
1 if kkF > 2.

The fermionic correction function WF (r) is defined in terms of gF (r) and

SF (k) as:
2 2 2√
~2 k 2 (SF (k) − 1) ~ ∇ gF

WF (k) = − (2SF (k) + 1) + F √ . (4.20)
4m SF (k) m gF
2 Now, the defect function is defined in Eq. (4.2). After taking the Fourier
transform of R(r), we find the new S(k):
1
S(k) = p , (4.21)
1 + 4γ/k 4 + (4/k 2 ) [(m/~2 )WF (k) − R(k)]
where γ = 4πn/a, and n is the electron density:

3 1
n= . (4.22)
4π (rs a)3
3 The new g(r) can be found directly as:
e−ik·r
Z
−1 1
g(r) − 1 = F [S(k) − 1] = d3 k [S(k) − 1]
n (2π)3
k 2 (S − 1) −1
Z Z Z 2π
4π 3 ikr cos θ
= − (rs a) dk d(cos θ)e dφ
3 (2π)3 1 0
2(rs a)3 3
Z
sin(kr)
≡ kF dk[S(k) − 1]k 2
3π kr
Z
3
= dk[S(k) − 1]k 2 sinc(kr).
2
(4.23)
In the third step of Eq. (4.23), we assigned k/kF as k.
38
Figure 4.6: The plot of the convergence of the radial distribution function
4 As explained in the 4th item of Sec. 4.1.2, we use mixing for each step and
iterate the system until it converges. Fig 4.6 exhibits the convergence of the
radial distribution function. To obtain this figure, we initiated the iteration
process with the function gF as defined in Eq. (4.18). To reach the conver-
gence of g(r), we did only 10 iterations, and used the same mixing factor
as Eq. (4.7). However, to calculate energy, we have done more iterations
to store the data for each rs value separately, and then interpolated them
to take the integral of the potential energy with respect to rs .
5 We find the potential energy for this case as:

Z
V n
= V (rs ) = d3 r(g(r) − 1)vc (r)
N 2
e2
Z
3 1 1 2
= r dr sin θdθdφ(g − 1) (4.24)
4π (r a)3 2 r
Zs
3
≡ rdr(g − 1) Ry,
rs
e2
where Ry= 2a
is the unit of the energy. The total energy per electron is:
Z rs
3 1
EGS = F + 2 drs0 rs0 V (rs0 ), (4.25)
5 rs 0
39
where F is the Fermi energy. It is found as the following:
2 2
~2 kF2 ~2 9π 3 1 ~2 9π 3 1 me2
F = = =
2m 2m 4 rs2 a2 2m 4 rs2 a ~2
" 2 # 2 (4.26)
e2 1 9π 3 1 9π 3
= = Ry.
2a 2 4 2 4
By placing the Fermi energy, the first term of the ground state energy is found
in terms of Ry as:
23
3 3 9π 2.2099
0 = F = ' (4.27)
5 10 4 rs2
The correlation energy of 3D electron gas:
Z rs
1
EC = −Ex + 2 drs0 rs0 V (rs0 ), (4.28)
rs 0
where Ex is the exchange energy per electron which is found as:
13
3 e2 kF

3 9π 0.916
Ex = − =− Ry ' − Ry. (4.29)
4 π 2πrs 4 rs
The unpolarized electron gas has only the Coulomb potential among the pair
of electrons, and obeys the Pauli exclusion principle. We are starting the itera-
tion process with the noninteracting case functions: gF (r), and its corresponding
SF (k). In the noninteracting case, electrons are so close to each other that we can
neglect the Coulomb interactions. This is the maximum density situation where
rs goes to zero. The pair distribution function can be defined as the probability
of finding another particle from our reference particle. Due to the Pauli principle,
the probability to find two electrons at the same place is 1/2. That is why in the
Fig. 4.7 we see that gF (r = 0) equals to 0.5 at the highest electron density level.
As the density decreases, one can see a decrease in the g(r) value at zero distance.
When the system reaches rs = 15, this value becomes almost zero, because it is
hard to find two electrons with opposite spins at the same point if they are not
in the high density regime.
40
Figure 4.7: Radial distribution function versus r/rs a for 3D electron gas
The correlation energy results of the 3D electron gas at different densities are
given in Table 4.3 along with the results of diffusion Monte Carlo, FHNC, and
coupled cluster method calculations. Although, DMC does not give ”exact” re-
sults for electrons, it is still a very reliable method. The difference between our
results and those of Ceperley et al. is typically less than 0.3% that we consider
fairly satisfactory for all practical purposes. In addition to this, as you can see
in the Table 4.3, we reproduced the results by Kallio et al., and confirmed that
this method works very well.
rs Ecpres EcKP EcDM C EcF HN C EcCCM EcV M C

1 -0.1155 -0.1159 -0.1220 -0.1141 -0.1220 -0.1855
2 -0.0869 -0.0870 -0.0902 -0.0859 -0.0904 -0.1050
3 -0.0718 -0.0719 -0.0738 -0.0710 -0.0738 -0.0807
4 -0.0621 -0.0620 -0.0636 -0.0612 -0.0634 -0.0675
5 -0.0550 -0.0549 -0.0563 -0.0541 -0.0560 -0.0585
10 -0.0362 -0.0360 -0.0372 -0.0355 -0.0370 -0.0372
20 -0.0225 -0.0224 -0.0230 -0.0218 -0.0236 -0.0228
Table 4.3: Correlation energy results for the unpolarized electron gas in Ry. KP
stands for the results of Kallio et al. [3]. DMC results from Ref. [18], FHNC
results from Ref. [19], CCM results from Ref.[20], and VMC results from Ref.[21].
The ground-state energies of Ref.[21] are converted to correlation energies for
comparison.
41
Since we have performed the computation from rs = 0 to rs = 1 by 0.05 steps, we
could find accurate results for high density regime, as you can see in the Table
4.4. If we increase the step numbers more, we can find a more accurate result
for rs = 0.1; because we are taking an integral from rs = 0 to rs = 0.1, and
increasing steps provides a better interpolation of the data. Table 4.4 shows that
for small rs , our numerical results agree closely with the field theory results [22].
Since field theoretical approach works at high electron density regime, we can
compare our small rs results with it.
rs Ecpres EcKP EcBL

0.1 -0.223 -0.229 -0.237
0.2 -0.190 -0.193 -0.194
0.3 -0.171 -0.172 -0.169
0.4 -0.157 -0.158 -0.151
0.5 -0.147 -0.147 -0.137
0.6 -0.138 -0.138 -0.126
0.7 -0.131 -0.132 -0.116
0.8 -0.125
0.9 -0.120
Table 4.4: Correlation energy results for the unpolarized electron gas in high
density range in Ry. The correlation energies ECBL are calculated from ECBL =
0.0622 ln rs − 0.094, Ref. [22], and ECKP are the results of Kallio et al. [3].
4.4 2D Spin Unpolarized Electron Liquid
4.4.1 Definition
The two-dimensional electron gas model can be realized in semi-conductor inter-

faces (GaAs-Ga1−x Alx As) or at the interface of a metal oxide and a semiconductor
(MOS) [23]. The electron liquid in such interfaces has important characteristics
such that they display the quantum Hall effect at very low temperatures and in
the presence of a strong and uniform magnetic field. Thereby, the theoretical
42
investigation of the behavior of 2D electron liquid has significance owing to its
intriguing, complex behavior and its several applications in many-body physics.
In this section, we present our results of the statistical properties for uniform 2D
electron gas at zero temperature with no magnetic field.
4.4.2 The Algorithm
In this part, the algoritm we used for 2D electron gas is explained step by step:
1 At first, we make all the variables dimensionless:
r/rs a → r, k/kF → k (4.30)

√
2
where kF is the Fermi wave vector. For 2D electron gas kF = rs a
. As
mentioned before, we write rs = 0 in the code for the starting case, and
write the known definitions of the noninteracting g(r) and S(k). For the
unpolarized 2D electron gas; these functions are:
2
1 2J1 (kF r)
gF (r) = 1 − . (4.31)
2 kF r
 q
 2 sin−1 k + k 1 − ( 2kkF )2 if k
≤2
π 2kF 2kF kF
SF (k) = (4.32)
1 k
if kF
> 2.
By using the definitions of gF (r) and SF (k), we can find the fermionic
correction term WF (r) as:
2 2 2√
~2 k 2 (SF (k) − 1) ~ ∇ gF

WF (k) = − (2SF (k) + 1) + F √ . (4.33)
4m SF (k) m gF
2 The defect function can be calculated by Eq. (4.2). After taking the Fourier
transform of R(r), we find the new S(k):
1
S(k) = p , (4.34)
1 + 4γ2d /k 3 + (4/k 2 ) [(m/~2 )WF (k) − R(k)]
43
where γ2d = 2πn/a, and n is the electron density:
1
n= . (4.35)
π(rs a)2
3 The new g(r) can be found directly as:
e−ik·r
Z
1
g(r) − 1 = d2 k [S(k) − 1]
n (2π)2
Z Z π
2 2 1
= πrs a dk(S − 1)k dθeikr cos θ (4.36)
(2π)2 0
Z
≡ dk[S(k) − 1]kJ0 (kr).
In the third step of the Eq. (4.36), we assigned k/kF as k.
4 We inserted new g(r) into the defect function R(r) which is defined in
Eq. (4.2). We did the iteration process until the pair distribution function
converges for each rs value. In the Fig. 4.8, you can see the convergence
of the radial distribution function. Each point on the plot stands for the
R
integral value dr|gaf ter (r) − gbef ore (r)| calculated after each iteration.
Figure 4.8: The plot of the convergence of the radial distribution function for 2D
electron gas
5 The potential energy for each rs value can be calculated by Eq. (4.16).
After taking interpolation of the function V (rs ), we can find the ground-
state energy:
44
Z rs
F 1
EGS = + 2 drs0 rs0 V (rs0 ), (4.37)
2 rs 0
where F is the Fermi energy. It is found as the following:
~2 kF2 ~2 2 ~2 2 me2
F = = =
2m 2m (rs a)2 2m (rs2 a) ~2
(4.38)
e2 2

2
= 2
= 2 Ry.
2a rs rs
By placing the Fermi energy, the first term of the ground state energy is found
in terms of Ry as:
F 1
0 = = 2. (4.39)
2 rs
The correlation energy of 2D electron gas:

Z rs
1
EC = −Ex + 2 drs0 rs0 V (rs0 ), (4.40)
rs 0
where Ex is the exchange energy per electron which is calculated as:
√
4 e2 kF 8 2 1.200
Ex = − =− Ry ' − Ry. (4.41)
3 π 3πrs rs
Our results of pair correlation function g(r) for 2D paramagnetic electron gas
for different rs are presented in Fig. 4.9. As we expected, the pair distribution
function at zero distance between the electrons approaches to zero as the density
of electrons decreases. In Sec. 4.2.3, we have stated that three-body correlations
are more important for 2D systems than 3D systems. However, by using the
simple and effective method by Kallio et al. for 2D, we could find good results
for the 2D unpolarized electron gas. To show this in a clear way, we compared
our present result g(r) with the FHNC result in the Fig. 4.10. Although the
result by Asgari et al. includes triplet correlations, it is obvious that there is not
a significant difference between them.
45
Figure 4.9: Radial distribution function of 2D electron gas at several densities for
unpolarized 2D electron gas.
Figure 4.10: The pair distribution function g(r) of the 2D electron gas at rs = 1 in
the paramagnetic state, as a function of distance r (in units of rs a). The present
results, shown in red, are compared with FHNC data by Asgari et al. [24], shown
in blue.
46
Bibliography
[1] G. F. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid.

Cambridge University Press, 2008.
[2] A. Fabrocini, S. Fantoni, and E. Krotscheck, Introduction to Modern Methods

of Quantum Many-Body Theory and Their Applications. World Scientific
Publishing, 2002.
[3] A. Kallio and J. Piilo, “Novel analytic calculation of electron gas properties,”
Physical Review Letters, vol. 77, no. 20, 1996.
[4] G. D. Mahan, Many Particle Physics. Physics of Solids and Liquids, 2002.
[5] J. M. J. V. Leeuwen, J. Groeneveld, and J. D. Boer, “New method for the

calculation of the pair correlation function. I,” Physica, vol. 25, pp. 663–667,
1959.
[6] S. Fantoni and S. Rosati, “Automatic classification for pathological prostate

images based on fractal analysis,” Nuovo Cimento A, vol. 20, no. 179, 1974.
[7] A. Kallio, V. Apaja, and S. Pöykkö, “Spectator fermion binding of bosons,”

Physica B, vol. 210, no. 472-478, 1995.
[8] V. Apaja, J. Halinen, V. Halonen, E. Krotscheck, and M. Saarela, “Charged-

boson fluid in two and three dimensions,” Physical Review B, vol. 55, no. 19,
1997.
[9] B. W. Ninham, “Charged Bose gas in astrophysics,” Physics Letters, vol. 4,

pp. 278–279, 1963.
47
[10] J. P. Hansen, B. Jancovici, and D. Schiff, “Phase diagram of a charged Bose
gas,” Physical Review Letters, vol. 29, no. 991, 1972.
[11] V. L. Ginzburg, “Superfluidity and superconductivity in the universe,” Jour-

nal of Statistical Physics, vol. 1, 1969.
[12] S. Moroni, S. Conti, and M. P. Tosi, “Monte carlo simulations of the charged
boson fluid at t=0,” Physical Review B, vol. 53, no. 9688, 1996.
[13] R. T. Hock-Kee Sim and F. Y. Wu, “Ground-state energy of charged quan-

tum fluids in two dimensions,” Physical Review B, vol. 34, no. 7123, 1986.
[14] A. Gold, “Local-field correction of the charged Bose condensate for two and
three dimensions,” Zeitschrift fr Physik B Condensed Matter, vol. 89, pp. 1–
10, 1992.
[15] M. Gell-Man and K. A. Brueckner, “Correlation energy of an electron gas

at high density,” Physical Review, vol. 106, no. 2, 1956.
[16] D. Ceperley, “Condensed-matter physics: Return of the itinerant electron,”

Nature, vol. 397, no. 386, 1999.
[17] N. H. March, Electron Correlations in the Solid State. Imperial College Press,
London, 1999.
[18] D. M. Ceperley and B. J. Alder, “Ground state of the electron gas by a

stochastic method,” Physical Review B, vol. 45, no. 566, 1980.
[19] L. J. Lantto, “Variational theory of multicomponent quantum fluids: An

application to positron-electron plasmas at t=0,” Physical Review B, vol. 36,
no. 5160, 1987.
[20] K. Emrich and J. G. Zabolitzky, “Electron correlations in the Bogoljubov

coupled-cluster formalism,” Physical Review B, vol. 30, no. 2049, 1984.
[21] D. Ceperley, “Ground state of the fermion one-component plasma: A Monte

Carlo study in two and three dimensions,” Physical Review B, vol. 18, no. 7,
pp. 3126–3138, 1978.
48
[22] R. F. Bishop and K. H. Lührmann, “Electron correlations: I. ground-state
results in the high-density regime,” Physical Review B, vol. 17, no. 3757,
1978.
[23] B.Tanatar and D. M. Ceperley, “Ground state of the two-dimensional elec-

tron gas,” Physical Review B, vol. 39, no. 8, 1988.
[24] R. Asgari, B. Davoudi, and M. Tosi, “Analytic theory of correlation energy

and spin polarization in the 2d electron gas,” Solid State Communications,
vol. 131, pp. 301–305, 2004.
49
Appendix A
Minimization of the Energy
As it is explained in Sec. 3.2, we will take functional derivative of the total energy
per particle expression (A.1) with respect to the function G(r). Since we want to
find the ground-state energy by a variational method, we equate the derivative
of the energy expression to zero in order to minimize the energy.
We have done the functional derivatives step by step to reach the iterative scheme
of the HNC/0 method. First, we write the definition of total energy per particle
in its open form:
~2 2 (S − 1)
3
Z Z
n 1
e= drg(r)v(r) − dkk
2 8m (2π)3 n S
2 Z (A.1)
n ~
− drg(r)∇2 ln g(r).
2 4m
p
As stated in Sec. 3.2 we take the functional derivative with respect to g(r)
which is denoted as G(r). First, we start from the first term of the energy
R
expression. The integral part of this term can be defined as F1 = drG2 (r)ν(r).
The functional derivative of an equation is defined by the following formula:
Z Z
δF ∂f ∂f
φ(r)dr = −∇ φ(r)dr, (A.2)
δn(r) ∂n ∂∇n
R
where F = drf (r). Here, F = F1 , and f (r) = G2 (r)v(r). So we can take the
50
functional derivative of the F1 function:
 zero

Z Z
∂ ∂
drG2 v(r) =  (G2 ν(r)) − ∇ G2 ν(r)  (δG)dr, (A.3)

δF1 = δ
∂G ∂∇G
Z
δF1
δF1 = 2G(r)ν(r)(δG)dr ⇒ = 2Gν(r). (A.4)
δG
Then the functional derivative of the first term is found as:
Z
δ n 2 n δF1 n
drG (r)ν(r) = = 2Gν(r) = nG(r)ν(r). (A.5)
δG 2 2 δG 2
Now, we can move our calculation to the second term of the energy per particle
expression, which is:
~2 (S − 1)3
Z
1
− dkk 2 , (A.6)
8m (2π)3 n S
where its integral part is denoted as F2 . Now, we take the functional derivative
of is as follows:
3
Z
2 (S − 1)
δF2 = δ dkk
S
Z 3 3

∂ 2 (S − 1) ∂ 2 (S − 1)
= k −∇ k (δS)dk (A.7)
∂S S ∂∇S S
∂ (S − 1)3
Z
= k2 δSdk.
∂S S
In this integral the partial derivative with respect to S(k) is found as:
2
∂ (S − 1)3

S−1
= (2S + 1). (A.8)
∂S S S
By placing Eq. (A.8) into (A.7) we can find;
(S − 1)2
Z
2
δF2 = k dk (2S + 1)δS, (A.9)
S
δ(S − 1) = δ(G2 − 1). (A.10)
In Eq. (A.9), δS can be found by using the Fourier transform relationship between
S − 1 and g − 1. Thus;
2−1
S − 1 = G^ (A.11)
2 − 1)
δ(S − 1) = δ (G^ (A.12)
51
]
δS = 2GδG (A.13)
By substituting Eq. (A.13) into the Eq. (A.9):

2
S−1
Z
2 ]
δF2 = k dk (2S + 1)2GδG (A.14)
S
Now, we can write the functional derivative of the second term of the energy
expression (A.6) by inserting Eq. (A.14) to it:
2 2
S−1 2 ]
Z
1 ~k
dk − (2S + 1)( ) GδG (A.15)
(2π)3 n 4m S
To solve Eq. (A.15), we will use Parseval’s identity:

Z Z
1 ] B(k),
]
n A(r)B(r)dr = dk A(k) (A.16)
(2π)3 n
] = − ~2 k2 (2S + 1)( S−1 )2 = w(k),

where A(k) ] and B(k)
] = GδG.
] Therefore, the
4m S
functional derivative of the second term with respect to G(r) can be found as:
~2 3

2 (S − 1)
Z
δ 1
− dkk , = nw0 (r)G(r) (A.17)
δG 8m (2π)3 n S
The integral part of the third term of the Eq. (A.1) is denoted as F3 :
Z Z
F3 [G] = dr(∇G) ⇒ F3 [G + h] = dr(∇G(r) + h0 (r))2
2
(A.18)
As → 0 the functional derivative of F3 :

Z
d d
δF3 = F[G + h] = dr(∇G + h0 )2 (A.19)
d d
The derivative with respect to is found as:
d d
(∇G + h0 )2 = (∇G)2 + 2∇Gh0 + 2 h02 = 2∇Gh0 (A.20)
d d
Then we put the Eq. (A.20) into (A.19):
Z
δF3 = dr2∇G∇h (A.21)
Now, we apply integration by parts:

Z Z Z
dr∇G∇h = G · h − h∇ Gdr = − h∇2 Gdr
2
(A.22)
52
So, δF3 is: Z
δF3 = −2 h∇2 Gdr. (A.23)
Then, we can write the functional derivative of F3 with respect to G(r0 ):

Z
δF3 [G]
= −2 ∇2 Gδ(r − r0 ) = −2∇2 G. (A.24)
δG(r0 )
Thus, we can place Eq. (A.24) into the third term, and then take its functional
derivative with respect to G:
n ~2 ~2 ~2
Z
δ
− drg(r)∆ ln g(r). = −2 n∇2 G = − n∇2 G(r). (A.25)
δG 2 4m 2m m
Now, we found the functional derivatives of all the terms of the energy per particle
expression given in Eq. (A.1). If place our results for 1st, 2nd and 3rd terms of
Eq. (A.1), and equate it to 0 for minimization of the energy, we find:
δe ~2
= nG(r)vc (r) + nw0 (r)G(r) − n∇2 G(r) = 0
δG 2 m
~ 2
n − ∇ G + vc (r) + w0 (r) G(r) = 0 (A.26)
m
2
~
− ∇2 G + vc (r) + w0 (r) G(r) = 0,
m
p
where G(r) = g(r). This is the main differential equation of the hypernetted
chain method. It is also called Schrödinger-like equation, because the function
G(r) acts like the probability amplitude of this zero energy equation. This is
explained in Sec. 3.4.
53
Appendix B
The Equivalence of the Methods
From the Eq. (3.22), the main differential equation is:
−~2 2
∇ ψ + (WB (r) + vc (r))ψ = 0, (B.1)
m
√
where ψ = g. Here, we will show the equivalence of the HNC/0 iterative schemes
which are presented by Kallio et al. and Fantoni. By this equivalence, we will
find how γ term is found as 4πn/a for 3D electron gas by Kallio, so that we can
convert it for 2D electron gas case.
The fluctuations of a function can cause crucial problems while taking its integral.
Since Kallio did the iteration on a defect function R(r), his method eliminates
this problem, and we can find better results. That’s why we need to find γ for
the 2D electron gas.
First we reorganize the highly nonlinear differential equation defined in Eq. (B.1):
−~2 ∇2 ψ
+ WB + vc = 0. (B.2)
m ψ
Then, we insert the definitions of vc (r) and WB (r) into the Eq. (B.2), and multiply
both sides by m/~2 :
54
" 2 #
∇2 ψ me2 k2 S − 1
+ 2 =F (2S + 1) . (B.3)
ψ ~r 4 S
Then, we take inverse Fourier transform of both sides:
2
k2 −∇2 ψ me2

S−1 −1
(2S + 1) = F + 2 . (B.4)
4 S ψ ~r
Now we put Eq. (3.27) into (B.4):
2
k2 me2

S−1 −1 1 2
(2S + 1) = F − ∇ (g − 1) − R(r) + 2 . (B.5)
4 S 2 ~r
Here, F −1 [∇2 (g−1)] = (−k 2 )(S−1) due to the Fourier transform relation between
g(r) − 1 and S(k) − 1. So, we insert this into Eq. (B.5). After multiplying both
sides by 4/k 2 , we find:
4me2

1 4 1
2S − 3 + 2 = (S − 1) − 2 R(k) + 2 2 F . (B.6)
S k ~k r
Here; ~2 /me2 = a is the Bohr radius and F[1/r] = 4πn/k 2 if we take the Fourier
transform in 3D analytically. Reorganizing the equation we find:

1 4 4 1 4 4 4πn
= 1 − 2 R(k) + 2 F = 1 − 2 R(k) + 2 , (B.7)
S2 k ak r k ak k2
−1/2
4 16πn
S(k) = 1 − 2 R(k) + . (B.8)
k ak 4
This is the main iterative formula which is also given in Eq. (3.28). From their
equivalence, one can find:
16πn 4γ
4
= 4. (B.9)
ak k
55
So, here we confirmed that γ = 4πn/a for 3D electron gas. To find 2D γ, we
return the Eq. (B.6) and define the Fourier transform for 2D system as:

1 2πn
F = (B.10)
r k
where n = 1/π(rs a)2 . By placing Eq. (B.10) into (B.6), we find the equation
1 4 4 2πn
2
= 1 − 2 R(k) + 2 (B.11)
S k ak k
So, the static structure factor is:
−1/2 −1/2
4 8πn 4 4γ2d
S= 1 − 2 R(k) + 3 = 1 − 2 R(k) + 3 (B.12)
k ak k k
where γ2d = 2πn/a.
56
Appendix C
The Code
In this chapter, I share the Mathematica code for 2D and 3D Bose and electron
liquid systems. The algorithms are explained in detail in Chapter 4. The choice
of range of integrals, interpolation, and the mixing ratio may need to be adapted
for different rs values and/or different cases.
C.1 3D Electron Gas
Below is the Mathematica code for the 3D spin unpolarized electron gas calcula-
tions.
(∗ Beginning o f c o n s t a n t s and f u n c t i o n d e f i n i t i o n s ∗)
i t e r c n t = 2 0 0 ; (∗ Number o f i t e r a t i o n s ∗)
r s = 1 ; (∗ Assign an i n i t i a l v a l u e f o r r s ∗)
g f [ r ] := 1 − ( ( ( 9 / 2 ) ∗ S p h e r i c a l B e s s e l J [ 1 , k f r [ r ] ] ˆ 2 ) / ( k f r [ r ] ˆ 2 ) ) ;
S f [ k ] := N[ I f [ k < 2 , (3∗ k ) /4 − ( k ˆ 3 ) / 1 6 , 1 ] ] ;
k f r [ r ] := ( ( ( 9 ∗ Pi ) / 4 ) ˆ ( 1 / 3 ) ) ∗ r ;
k o k g f [ r ] := Sqrt [ g f [ r ] ]
kokgr [ r ] := Sqrt [ g r [ r ] ]
Lapgminus [ r ] := Derivative [ 2 ] [ g r ] [ r ] + (2∗ Derivative [ 1 ] [ g r ] [ r ] ) / r ;
57
Lapkokgf [ r ] :=
Derivative [ 2 ] [ k o k g f ] [ r ] + (2∗ Derivative [ 1 ] [ k o k g f ] [ r ] ) / r ;
Lapkokgr [ r ] :=
Derivative [ 2 ] [ kokgr ] [ r ] + (2∗ Derivative [ 1 ] [ kokgr ] [ r ] ) / r ;
ad [ r ] := ( Lapkokgf [ r ] / k o k g f [ r ] ) − ( Lapkokgr [ r ] / kokgr [ r ] ) + ( 1 / 2 ) ∗
Lapgminus [ r ] ;
adFT0 [ r , k ] := ( r ˆ2)∗ ad [ r ] ∗ S i n c [ k∗ k f r [ r ] ] ;
add [ k ] := 3∗
NIntegrate [ adFT0 [ r , k ] , { r , 0 , 4 9 . 9 1 } ,
Method −> {” GlobalAdaptive ” , ” S y m b o l i c P r o c e s s i n g ” −> 0 } ] ;
Smin [ k ] := ( ( 1 + ( ( 4 / ( 9 ∗ Pi ) ) ˆ ( 4 / 3 ) ) ∗ ( 1 2 ∗
r s /( k ˆ4)) + ((2∗ Sf [ k ] +
1 ) ∗ ( ( ( S f [ k ] − 1 )/
S f [ k ] ) ˆ 2 ) ) + ( ( ( ( 4 ∗ ( ( 4 / ( 9 ∗ Pi ) ) ˆ ( 2 / 3 ) ) ) ) / ( k ˆ 2 ) ) ∗
add [ k ] ) ) ˆ ( − 1 / 2 ) ) − 1 ;
SminFT0 [ k , r ] := ( k ˆ2)∗ S i n t e r [ k ] ∗ S i n c [ k∗ k f r [ r ] ] ;
smin [ r ] := ( 3 / 2 ) ∗
NIntegrate [ SminFT0 [ k , r ] , {k , 0 , 4 9 . 9 1 } ,
Method −> {” GlobalAdaptive ” , ” S y m b o l i c P r o c e s s i n g ” −> 0 } ]
(∗ End o f d e f i n i t i o n s ∗)
(∗ The b l o c k b e l o w r e p r e s e n t s a s i n g l e i t e r a t i o n . I t t a k e s a t r i a l
f u n c t i o n as an i n p u t , c a l c u l a t e s S ( k ) and u s i n g t h a t S ( k ) ( s i n t e r ) ,
c a l c u l a t e d g ( r ) ( g i n t e r ) . e r r i s t h e v a r i a t i o n between t h e t r i a l
f u n c t i o n and t h e c a l c u l a t e d f u n c t i o n . Output f u n c t i o n g o u t i s t h e
m i x t u r e o f g t r i a l and g i n t e r ∗)
E r r o r C a l c g r [ g r t r i a l ] := Block [ { } ,
gr = g r t r i a l ;
S i n t e r = Quiet [ FunctionInterpolation [ Smin [ k ] , {k , 0 . 0 1 , 4 9 . 9 1 } ] ] ;
ginter =
Quiet [ FunctionInterpolation [ ( 1 + smin [ r ] ) , { r , 0 . 0 1 , 4 9 . 9 1 } ] ] ;
(∗ t o r e d u c e t h e f l u c t u a t i o n s c a l c u l a t e d g ( r ) i s mixed w i t h t h e
t r i a l f u n c t i o n ∗)
58
f i f [ r ] := ( ( g i n t e r [ r ] ) + (9∗ g r t r i a l [ r ] ) ) / 1 0 ;
gout = Interpolation [
Transpose [ { Table [ n , {n , 0 . 0 1 , 4 9 . 9 1 , 0 . 2 } ] ,
Table [ f i f [ n ] , {n , 0 . 0 1 , 4 9 . 9 1 , 0 . 2 } ] } ] ] ;
e r r = NIntegrate [ Abs [ gout [ r ] − g r t r i a l [ r ] ] , { r , 0 , 1 2 } ] ;
Return [ { e r r , gr , gout , g i n t e r } ] ;
]
(∗ The b l o c k b e l o w do es t h e i t e r a t i o n s f o r a g i v e n r s and c o l l e c t s t h e
r e s u l t o f each i t e r a t i o n ∗)
Rstogr [ r s n ] := Block [ { } ,
i t e r = {};
Do[
rs = rsn ;
e r r 1 = ErrorCalcgr [ gr ] ;
i t e r = Append [ i t e r , e r r 1 ] ;
Print [ ToString [ n ] <> ” E r r o r : ” <> ToString [ Part [ i t e r , n , 1 ] ] ] ;
g r = Part [ i t e r , n , 3 ] ;
I f [ e r r < 0 . 0 0 0 0 0 1 , Break [ Null , Do ] ]
, {n , i t e r c n t } ] ; i t e r ] ;
C.2 3D Charged Bose Gas
Below is the Mathematica code for the 3D Charged boson gas calculations.
i t e r c n t = 250;
rs = 0;
g f [ r ] := 1 ;
S f [ k ] := 1 ;
k f r [ r ] := ( ( ( 9 ∗ Pi ) / 4 ) ˆ ( 1 / 3 ) ) ∗ r ;
kokgr [ r ] := Sqrt [ g r [ r ] ] ;
Lapgminus [ r ] := Derivative [ 2 ] [ g r ] [ r ] + (2∗ Derivative [ 1 ] [ g r ] [ r ] ) / r ;
Lapkokgr [ r ] :=
59
Derivative [ 2 ] [ kokgr ] [ r ] + (2∗ Derivative [ 1 ] [ kokgr ] [ r ] ) / r ;
ad [ r ] := ( Lapkokgr [ r ] / kokgr [ r ] ) − ( 1 / 2 ) ∗ Lapgminus [ r ] ;
adFT0 [ r , k ] := ( r ˆ2)∗ ad [ r ] ∗ S i n c [ k∗ k f r [ r ] ] ;
add [ k ] := 3∗
NIntegrate [ adFT0 [ r , k ] , { r , 0 , 10 } ,
Smin [ k ] := ( ( 1 + ( ( 4 / ( 9 ∗ Pi ) ) ˆ ( 4 / 3 ) ) ∗ ( 1 2 ∗
r s / ( k ˆ 4 ) ) − ( ( ( ( 4 ∗ ( ( 4 / ( 9 ∗ Pi ) ) ˆ ( 2 / 3 ) ) ) ) / ( k ˆ 2 ) ) ∗ add [ k ] ) ) ˆ ( − 1 /
2)) − 1;
SminFT0 [ k , r ] := ( k ˆ2)∗ S i n t e r [ k ] ∗ S i n c [ k∗ k f r [ r ] ] ;
gmin [ r ] := ( 3 / 2 ) ∗
NIntegrate [ SminFT0 [ k , r ] , {k , 0 , 40} ,
gr = gf ;
S i n t e r = Quiet [ FunctionInterpolation [ Smin [ k ] , {k , 0 , 4 0 } ] ] ;
g i n t e r = Quiet [ FunctionInterpolation [ ( 1 + gmin [ r ] ) , { r , 0 , 4 0 } ] ] ;
f i f [ r ] := ( ( g i n t e r [ r ] ) + (2∗ g r t r i a l [ r ] ) ) / 3 ;
Transpose [ { Table [ n , {n , 0 . 0 1 , 4 0 , 0 . 1 } ] ,
Table [ f i f [ n ] , {n , 0 . 0 1 , 4 0 , 0 . 1 } ] } ] ] ;
e r r = NIntegrate [ Abs [ gout [ r ] − g r t r i a l [ r ] ] , { r , 0 , 1 2 } ] ;
]
i t e r = {};
Do[
rs = rsn ;
60
g r = Part [ i t e r , n , 3 ] ;
I f [ e r r < 0 . 0 0 0 0 0 0 1 , Break [ Null , Do ] ]
, {n , i t e r c n t } ] ;
iter ];
C.3 2D Charged Bose Gas
Below is the Mathematica code for the 2D Charged boson gas calculations.
i t e r c n t = 150;
rs = 0;
g f [ r ] := 1 ;
gr = gf ;
grminus [ r ] := g r [ r ] − 1 ;
S f [ k ] := 1
k f r [ r ] := Sqrt [ 2 ] ∗ r ;
kokgr [ r ] := Sqrt [ g r [ r ] ] ;
Lapgminus [ r ] :=
Derivative [ 2 ] [ grminus ] [ r ] + ( Derivative [ 1 ] [ grminus ] [ r ] ) / r ;
Lapkokgr [ r ] := Derivative [ 2 ] [ kokgr ] [ r ] + ( Derivative [ 1 ] [ kokgr ] [ r ] ) / r ;
R[ r ] := ( Lapkokgr [ r ] / kokgr [ r ] ) − ( 1 / 2 ) ∗ Lapgminus [ r ] ;
RFT0 [ r , k ] := r ∗R[ r ] ∗ S i n c [ k∗ k f r [ r ] ] ;
Rk [ k ] :=
2∗ NIntegrate [ RFT0 [ r , k ] , { r , 0 , 9 9 . 9 1 } ,
Smin [ k ] := ( ( 1 + ( Sqrt [ 8 ] ∗ r s / ( k ˆ 3 ) ) − ( ( 2 / ( k ˆ 2 ) ) ∗Rk [ k ] ) ) ˆ ( − 1 / 2 ) ) − 1 ;
SminFT0 [ k , r ] := k∗ S i n t e r [ k ] ∗ S i n c [ k∗ k f r [ r ] ] ;
gmin [ r ] :=
NIntegrate [ SminFT0 [ k , r ] , {k , 0 , 9 9 . 9 1 } ,
61
Sinter =
Interpolation [
Table [ Smin [ n ] , {n , 0 . 0 1 , 9 9 . 9 1 , 0 . 1 } ] } ] ] ;
ginter =
Interpolation [
Table [ 1 + gmin [ n ] , {n , 0 . 0 1 , 9 9 . 9 1 , 0 . 1 } ] } ] ] ;
f i f [ r ] := ( ( g i n t e r [ r ] ) + (19∗ g r t r i a l [ r ] ) ) / 2 0 ;
Table [ f i f [ n ] , {n , 0 . 0 1 , 9 9 . 9 1 , 0 . 1 } ] } ] ] ;
e r r = NIntegrate [ Abs [ g i n t e r [ r ] − g r t r i a l [ r ] ] , { r , 0 , 9 9 . 9 1 } ] ;
]
i t e r = {};
Do[
rs = rsn ;
g r = Part [ i t e r , n , 3 ] ;
I f [ e r r < 0 . 0 0 0 0 0 0 3 8 6 , Break [ Null , Do ] ]
, {n , i t e r c n t } ] ;
iter ];
62
C.4 2D Electron Gas
Below is the Mathematica code for the 2D spin unpolarized electron gas calcula-
tions.
i t e r c n t = 120;
rs = 0;
g f [ r ] := 1 − ( ( 1 / 2 ) ∗ ( 2 ∗ BesselJ [ 1 , k f r [ r ] ] / k f r [ r ] ) ˆ 2 ) ;
S f [ k ] := N[
I f [ k <= 2 , ( 2/ Pi ) ∗ ( ArcSin [ k / 2 ] + ( k /2)∗ Sqrt [ 1 − ( ( k / 2 ) ˆ 2 ) ] ) , 1 ] ] ;
k f r [ r ] := Sqrt [ 2 ] ∗ r ;
k o k g f [ r ] := Sqrt [ g f [ r ] ]
kokgr [ r ] := Sqrt [ g r [ r ] ]
Lapgminus [ r ] := Derivative [ 2 ] [ g r ] [ r ] + ( Derivative [ 1 ] [ g r ] [ r ] ) / r ;
Lapkokgf [ r ] := Derivative [ 2 ] [ k o k g f ] [ r ] + ( Derivative [ 1 ] [ k o k g f ] [ r ] ) / r ;
Lapkokgr [ r ] := Derivative [ 2 ] [ kokgr ] [ r ] + ( Derivative [ 1 ] [ kokgr ] [ r ] ) / r ;
ad [ r ] := ( Lapkokgf [ r ] / k o k g f [ r ] ) − ( Lapkokgr [ r ] / kokgr [ r ] ) + ( 1 / 2 ) ∗
Lapgminus [ r ] ;
adFT0 [ r , k ] := r ∗ad [ r ] ∗ BesselJ [ 0 , k∗ k f r [ r ] ] ;
add [ k ] := 2∗
Quiet [ NIntegrate [ adFT0 [ r , k ] , { r , 0 , 1 5 . 0 1 } ,
Method −> {” GlobalAdaptive ” , ” S y m b o l i c P r o c e s s i n g ” −> 0 } ] ] ;
Smin [ k ] := ( ( 1 + ( ( Sqrt [ 8 ] ) ∗
r s /( k ˆ3)) + ((2∗ Sf [ k ] +
1 ) ∗ ( ( ( S f [ k ] − 1 )/ S f [ k ] ) ˆ 2 ) ) + ( 2 / ( k ˆ 2 ) ) ∗ add [ k ] ) ˆ ( − 1 / 2 ) ) − 1 ;
SminFT0 [ k , r ] := k∗ S i n t e r [ k ] ∗ BesselJ [ 0 , k∗ k f r [ r ] ] ;
Gmin [ r ] :=
Quiet [ NIntegrate [ SminFT0 [ k , r ] , {k , 0 , 1 5 . 0 1 } ,
Method −> {” GlobalAdaptive ” , ” S y m b o l i c P r o c e s s i n g ” −> 0 } ] ]
gr = gf ;
Sinter =
63
Interpolation [
Table [ Smin [ n ] , {n , 0 . 0 1 , 1 5 . 0 1 , 0 . 1 } ] } ] ] ;
ginter =
Interpolation [
Table [ 1 + Gmin [ n ] , {n , 0 . 0 1 , 1 5 . 0 1 , 0 . 1 } ] } ] ] ;
e r r = NIntegrate [ Abs [ g i n t e r [ r ] − g r t r i a l [ r ] ] , { r , 0 , 1 0 } ] ;
f i f [ r ] := ( g i n t e r [ r ] + ( 9 ∗ ( g r t r i a l [ r ] ) ) ) / 1 0 ;
Table [ f i f [ n ] , {n , 0 . 0 1 , 1 5 . 0 1 , 0 . 1 } ] } ] ] ;
]
i t e r = {};
Do[
rs = rsn ;
g r = Part [ i t e r , n , 3 ] ;
I f [ Part [ i t e r , n , 1 ] < 0 . 0 0 0 0 0 4 9 6 2 7 3 , Break [ Null , Do ] ]
, {n , i t e r c n t } ] ;
iter ];
64

Analytic Calculation of Ground State Properties of The 2D and 3D Electron Gas

Uploaded by

Copyright:

Available Formats

Analytic Calculation of Ground State Properties of The 2D and 3D Electron Gas

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Analytic Calculation of Ground State Properties of The 2D and 3D Electron Gas

Uploaded by

Copyright:

Available Formats

ANALYTIC CALCULATION OF GROUND

STATE PROPERTIES OF THE 2D AND 3D

Asst. Prof. Dr. Balázs Hetényi(Advisor)

Prof. Dr. Bilal Tanatar

Assoc. Prof. Dr. Hande Toffoli

Approved for the Graduate School of Engineering and Science:

Prof. Dr. Levent Onural

2.1 Jellium Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1.1 Particle Density and Wigner-Seitz Radius . . . . . . . . . 4

2.2 Coulomb Interaction . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.3 Radial Distribution Function . . . . . . . . . . . . . . . . . . . . . 6

2.3.1 Static Structure Factor . . . . . . . . . . . . . . . . . . . . 8

2.4 Description of a Quantum Many Body System . . . . . . . . . . . 9

3 Hypernetted Chain Method 13

3.2 Minimization of the Energy . . . . . . . . . . . . . . . . . . . . . 20

3.3 Variational Approach for Fermi Systems . . . . . . . . . . . . . . 22

3.4 HNC/0 Approximation For Fermions . . . . . . . . . . . . . . . . 23

4 Two and Three Dimensional Many Body Systems 26

4.1 3D Charged Boson Gas . . . . . . . . . . . . . . . . . . . . . . . . 26

4.1.2 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 27

4.1.3 Results and Analysis . . . . . . . . . . . . . . . . . . . . . 30

4.2 2D Charged Boson Fluid . . . . . . . . . . . . . . . . . . . . . . . 32

4.2.2 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 33

4.2.3 Results and Analysis . . . . . . . . . . . . . . . . . . . . . 34

4.3 3D Spin Unpolarized Electron Gas . . . . . . . . . . . . . . . . . 36

4.3.2 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 37

4.3.3 Results and Analysis . . . . . . . . . . . . . . . . . . . . . 40

4.4 2D Spin Unpolarized Electron Liquid . . . . . . . . . . . . . . . . 42

4.4.2 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 43

4.4.3 Results and Analysis . . . . . . . . . . . . . . . . . . . . . 45

A Minimization of the Energy 50

B The Equivalence of the Methods 54

C.1 3D Electron Gas . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

C.2 3D Charged Bose Gas . . . . . . . . . . . . . . . . . . . . . . . . 59

C.3 2D Charged Bose Gas . . . . . . . . . . . . . . . . . . . . . . . . 61

C.4 2D Electron Gas . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

3.1 The representation of cluster functions . . . . . . . . . . . . . . . 15

3.2 The notations of classified cluster diagrams . . . . . . . . . . . . . 17

3.3 The diagrammatic representation of the operation between nodal

3.4 An elementary diagram sample with four point structure. . . . . . 20

4.1 The plot of the convergence of the radial distribution function . . 29

4.2 Radial distribution function versus r/rs a for 3D Boson gas . . . . 30

4.3 Static structure factor versus kr0 for 3D Bose gas . . . . . . . . . 31

4.5 Radial distribution function versus r/rs a for 2D Boson gas . . . . 35

4.6 The plot of the convergence of the radial distribution function . . 39

4.7 Radial distribution function of 3D electron gas . . . . . . . . . . . 41

4.9 Radial distribution function of 2D electron gas at several densities

4.10 Radial distribution function of 2D electron gas compared with

4.1 Ground-state energy of the 3D charged boson fluid in Ry . . . . . 32

4.2 Ground-state energy results for 2D charged boson gas in Ry . . . 36

4.3 Correlation energy results for the unpolarized electron gas in Ry . 41

After science has stepped to quantum mechanics from classical mechanics,

In Chapter 3 we present the HNC Method aimed to investigate the distribution of

Our approximations are:

2.1 Jellium Model

2.1.1 Particle Density and Wigner-Seitz Radius

2.2 Coulomb Interaction

We can define the one-electron functions as plane waves:

g(r) = g (r) + g↑↓ (r). (2.10)

where F is the Fermi energy. It is found as the following: