VARIATIONAL MONTE CARLO
CALCULATIONS FOR BOSE-HUBBARD
MODEL BASED ON PROJECTED
WAVEFUNCTIONS
a thesis
submitted to the department of physics
and the graduate school of engineering and science
of bilkent university
in partial fulfillment of the requirements
for the degree of
master of science
By
Fulya Ko¸c
June, 2014
I certify that I have read this thesis and that in my opinion it is fully adequate, in scope and in quality, as a thesis for the degree of Master of Science.
Asst. Prof. Dr. Balazs Het´enyi(Advisor)
I certify that I have read this thesis and that in my opinion it is fully adequate, in scope and in quality, as a thesis for the degree of Master of Science.
Prof. Dr. Cemal Yalabık
I certify that I have read this thesis and that in my opinion it is fully adequate, in scope and in quality, as a thesis for the degree of Master of Science.
Assoc. Prof. Dr. Azer Kerimov
Approved for the Graduate School of Engineering and Science:
Prof. Dr. Levent Onural Director of the Graduate School
ABSTRACT
VARIATIONAL MONTE CARLO CALCULATIONS
FOR BOSE-HUBBARD MODEL BASED ON
PROJECTED WAVEFUNCTIONS
Fulya Ko¸cM.S. in Physics
Supervisor: Asst. Prof. Dr. Balazs Het´enyi June, 2014
Bose-Hubbard model is mainly used to describe and study the interactions be-tween neutral atomic gases trapped in an optical lattice [1] and Josephson junction arrays [2]. It is one of the toy models to understand quantum phase transitions, i.e. a phase transition exists between the Mott insulator state and the super-fluid state. Analytical solutions are limited to obtaining the ground state energy for small systems, whereas, computational studies can be done for larger system sizes. We applied the variational Monte Carlo method to the Bose-Hubbard model based on projected wavefunctions, i.e. Baeriswyl and Gutzwiller-Baeriswyl. Even though our method can be applicable to any dimension, we only consider the one dimensional case in this thesis. We expressed observables in forms of averages over configurations to which we can apply Monte Carlo sampling techniques. Our results for both Baeriswyl and Gutzwiller projections are in qualitatively good agreement with the known calculations of the phase diagram [3, 4]. Furthermore, we introduced a new method, apart from other known methods [5, 6], based on the Drude weight [7–9] to calculate the superfluid fraction, which can also be extended to observe BCS superconductivity [10].
¨
OZET
BOSE-HUBBARD MODEL˙I ˙IC
¸ ˙IN ˙IZD ¨
US
¸ ¨
UML ¨
U DALGA
FONKS˙IYONLARI KULLANARAK VARYASYONEL
MONTE CARLO HESAPLAMALARI
Fulya Ko¸cFizik, Y¨uksek Lisans
Tez Y¨oneticisi: Asst. Prof. Dr. Balazs Het´enyi Haziran, 2014
Bose-Hubbard modeli, genel anlamda, optik kafes i¸cerisine sıkı¸smı¸s n¨otr gazlar arasındaki etkile¸simler [1] ile Josephson eklemlerini [2] incelemek i¸cin kullanılır. Bose-Hubbard Hamiltonian’ı, kuvantum faz ge¸ci¸slerini, Mott yalıtkanı fazı ile s¨uperakı¸skan fazı arası ge¸ci¸si, anlamak i¸cin en temel modeli olu¸sturur. Anali-tik ¸c¨oz¨umler k¨u¸c¨uk sistemlerle sınırlı kalırken hesaplamalı ¸c¨oz¨umler daha b¨uy¨uk sistemleri anlamak i¸cin kullanılır. Bu tezde Bose-Hubbard modeli izd¨u¸s¨umsel dalga fonksiyonları, Baeriswyl ve Gutzwiller-Baeriswyl, kullanarak varyasyonel Monte Carlo yakla¸sımı ile incelenmi¸stir. Geli¸stirdi˘gimiz metod herhangi bir boyuta uygulanabilmesine ra˘gmen, bu tezde, sadece 1 boyutta hesaplamalar yapılmı¸stır. Fiziksel nicelikler Monte Carlo tekniklerini uygulayabilece˘gimiz ¸sekilde konfig¨urasyonlar ¨uzerinden ortalama hesapları yapılarak incelenmi¸stir. Elde etti˘gimiz sonu¸clar, daha ¨onceki bulgularla [3,4] nitelik bakımından iyi ¸sekilde ¨ort¨u¸smektedir. Ayrıca s¨uperakı¸skan fazı oranını hesaplamak i¸cin daha ¨once kul-lanılan metodlardan [5, 6] daha farklı bir metod, Drude a˘gırlı˘gı [7–9], i¸cin temel olu¸sturduk. Aynı zamanda bu metod BCS s¨uperiletkenli˘gini a¸cıklamak i¸cin de geni¸sletilebilir [10].
Anahtar s¨ozc¨ukler : Bose-Hubbard modeli, Gutzwiller, Baeriswyl, Varyasyonel
Contents
1 Introduction 1
1.1 Bose-Hubbard Model . . . 1
1.2 Symmetries of The Bose-Hubbard Hamiltonian . . . 3
1.3 Limiting Cases of the Bose-Hubbard Hamiltonian . . . 4
1.3.1 Mott Insulating Phase . . . 4
1.3.2 Superfluid Phase . . . 5
1.4 Outline of the thesis . . . 6
2 Background and General Definitions 7 2.1 Bose-Einstein Condensation . . . 7
2.2 Superfluid Phase . . . 10
3 Bose Hubbard Methods 18 3.1 Analytical Approach . . . 18
3.1.1 Mean-Field Theory . . . 18
CONTENTS vi
3.2 Computational Approach . . . 24
3.2.1 Exact Diagonalisation Method . . . 24
3.2.2 Quantum Monte Carlo Method . . . 26
3.2.3 Density Matrix Renormalization Group Method . . . 28
4 Algorithm 32 4.1 Variational Monte Carlo Method . . . 32
4.1.1 Metropolis-Hastings Algorithm . . . 33
4.2 Adaptation of the VMC Method to the Bose-Hubbard Model . . . 34
4.2.1 Baeriswyl Projection . . . 34
4.2.2 Gutzwiller-Baeriswyl Projection . . . 34
4.2.3 Calculation of The Observables . . . 35
5 Data Analysis and Conclusion 41 5.1 Baeriswyl Projection Results . . . 42
5.2 Gutzwiller-Baeriswyl Projection Results . . . 46
List of Figures
5.1 Propagator vs lattice site for Baeriswyl with L = 40 . . . 43
5.2 ET OT vs α for Baeriswyl with L = 40 . . . 44
5.3 Phase diagram for Baeriswyl with L=40 . . . 45
5.4 Phase diagram for Baeriswyl with L=20 . . . 46
5.5 ET OT vs α for Gutzwiller-Baeriswyl with L = 40 . . . 47
5.6 Phase diagram for Gutzwiller-Baeriswyl with L=40 . . . 48 5.7 Phase diagram comparison between GW and GW-BW with L = 40 49
Chapter 1
Introduction
1.1
Bose-Hubbard Model
The Bose-Hubbard model is a bosonic analogue of the Fermi-Hubbard model [11] which studies the strongly correlated materials. The Bose-Hubbard model was first introduced by Gersch and Knollman in 1963 [12], and further studied exten-sively by Fisher et al. in 1989 [13]. Based on the analysis of Fisher et al., the phase diagram of the Bose-Hubbard model at zero temperature contains two dif-ferent phases, namely; Mott insulating phase and superfluid phase. Hamiltonian of the model is ˆ H = −t X <i,j> ˆb† iˆbj + h.c. + U 2 X i ˆ ni(ˆni− 1) − µ X i ˆ ni (1.1)
where ˆb†i and ˆbi boson creation and annihilation operators on site-i respectively.
ˆ
ni = ˆb†iˆbi is boson occupation number at site-i, and t is the hopping parameter.
< i, j > indicates that the summation over nearest neighboring sites. U, denotes on-site interaction; it can be either repulsive, U > 0, or attractive, U < 0. Finally, µ is the chemical potential.
Boson creation and annihilation operators obey the following commutation rela-tions; h ˆbi, ˆb†j i = δij and h ˆb† i, ˆb † i i =hˆbi, ˆbi i = 0 ∀i, j (1.2)
These operators also act on the eigenbasis of the single-site boson occupation number operator as;
ˆb†
i|ˆni > =
√
ni+ 1|ni+ 1 >
ˆbi|ˆni > = √ni|ni− 1 > (1.3)
The first term in (1.1) denotes the hopping term and it introduces hoppings of bosons between nearest neighboring sites. This term of the Hamiltonian, in a way, describes how particles are delocalized; hence, it is convenient to represent it as the kinetic term as well. Second term in (1.1), on the other hand, is called on-site
repulsion term, which contrarily to the first term tries to localize the bosons on
each site. The last term in (1.1) controls the particle number in the system. In order to analyse the quantum phase transition and the existence of two phases as indicated at the beginning of this chapter, we need to consider two cases; first one is Ut → 0 which corresponds to the localized case, i.e. Mott insulating phase, and the second one is t
U → ∞ which corresponds to a delocalized case, i.e.
superfluid phase. However before analysing these two limits, one needs to define the symmetries of the Bose-Hubbard Hamiltonian.
1.2
Symmetries of The Bose-Hubbard
Hamilto-nian
Basically, we can analyse the symmetry consideration of the Bose-Hubbard model in three different ways [14–16];
1. U (1) Symmetry:
This symmetry indicates the conservation of total number of particles in the system, and ˆH is invariant under transformation;
ˆb† i, ˆbi →ˆb†ie iθ, ˆb ie−iθ = ei ˆN θˆb† i, ˆbi e−i ˆN θ ∀θ ∈ ℜ (1.4) 2. Translational Symmetry:
This symmetry indicates the conservation of total quasi-momentum of par-ticles in the system, and ˆH is invariant under transformation;
ˆb† i+1, ˆbi+1 → e−i ˆT ˆb†i, ˆbi ei ˆT (1.5) with ˆT being; ˆ T = L−1 X k=0 2πk L ˆb † iˆbi (1.6)
where L is the lattice site. 3. Reflection Symmetry:
In this symmetry ˆH is invariant under transformation; ˆb† i, ˆbi →ˆb†N−i, ˆbN−i in real space or ˆb† k, ˆbk →ˆb†−k, ˆb−k in momentum space. (1.7)
1.3
Limiting Cases of the Bose-Hubbard
Hamil-tonian
1.3.1
Mott Insulating Phase
This limit corresponds to t
U → 0, and the Bose-Hubbard Hamiltonian reduces to;
ˆ Hon−site= U 2 X i ˆ ni(ˆni− 1) − µ X i ˆ ni (1.8)
Hence, the ground state wavefunction is just the product state of single particles, which can be represented as [17];
|ΨM I >= Ld Y i ˆb† i n |0 > (1.9)
where n is the boson per site, L is lattice site, and d is the dimensionality. (1.8) is minimized for integer values of n0 = 12 +Uµ. As we have a commensurate
fillings for a finite range of the chemical potential, we can say that the ground state, in this case, is incompressible, where compressibility is defined as κ = ∂ρ∂µ with ρ is the boson density [17].
As we consider the particle correlations in the ground state, we can realize that expectation value for such a correlation in momentum space, for practical reasons, i.e. < ˆb†
qˆbq > is independent of the value q; meaning, we have a delocalization in
momentum space as each momentum has the same weight; contrarily, we have localization in real space.
1.3.2
Superfluid Phase
This limit corresponds to Ut → ∞, and the Bose-Hubbard Hamiltonian reduces to; ˆ Hhopping = −t X <i,j> ˆb† iˆbj+ h.c. (1.10)
Since, hopping term is diagonal in momentum space, we can further represent it as; ˆ Hhopping = X ki ˆ nkiǫ(ki) with ǫ(ki) = −2t cos (kixi) (1.11)
For such a Hamiltonian, we cannot use the same ground state wavefunction any-more; however, we can define a new one as [18];
|ΨN >U=0 = 1 √ N ! 1 √ L! X − →x ˆb† − →x !N |0 > (1.12)
Since all particles would occupy the possible lowest energy, ground state of the Hamiltonian would be at k = 0; meaning, particles are well localized in a single state with a coherent phase [17] in momentum space; whereas, they are delocalized in real space.
Observe that since the bosons, now, have coherent phase, particle number is allowed to fluctuate, which breaks the U (1) symmetry as mentioned in Subsec-tion.1.2. This broken-symmetry state is called superfluid state.
1.4
Outline of the thesis
In this thesis, our goal is to analyse one dimensional Bose-Hubbard model by using projected wavefunctions; i.e. namely Baeriswyl, and Gutzwiller-Baeriswyl, via variational Monte Carlo calculations. We examine the observables in forms of averages over configurations which we can apply MC techniques. Critical values of the Mott lobes existing in the phase diagrams for the two different variational wavefunctions, which we used in this thesis, are analysed and compared with each other, and also with a reference value [19]. Besides, we calculated the superfluid fraction with a new method [9,10] taking Drude weight as basis, apart from other known methods [5, 6].
The thesis is organized as; in Chapter 2, some basic background and general definitions which are deeply related to the Bose-Hubbard model are given, in
Chapter 3, both analytical and computational methods done so far are analysed
basically, in Chapter 4, the algorithm we developed is introduced, and finally in the last chapter, results that we obtained from the projected wavefunctions are presented and compared.
Chapter 2
Background and General
Definitions
2.1
Bose-Einstein Condensation
Consider a many-body wavefunction which is symmetric under exchange of pairs i and j; Ψ (~r1, ..., ~rN) where N is the number of bosons. The single-particle density
matrix can be written as [20];
ρ1 ~r, ~r′ ≡ NX i pi Z d~r2...d ~rNΨ∗i(~r, ~r2, ..., ~rN)Ψi(~r′, ~r2, ..., ~rN) ≡ < ˆΨ†(~r) ˆΨ(~r′ ) > (2.1)
where pi is the probability of state i. Density matrix can be further expanded to
include spin and time as well.
Since ρ1 matrix is hermitian, i.e. ρ1(~r, ~r′) ≡ ρ∗1(~r
′
, ~r), we can diagonalise it in such a form [20];
ρ1 ~r, ~r′ =X i niχ∗i(~r)χi(~r′) (2.2)
where χi(~r) forms an orthonormal set, and nis are being the eigenvalues of this
set.
Based on (2.2), we can give a formal definition of BEC as [20];
1. BEC does not occur if all the eigenvalues of (2.2) are of the order unity. 2. Simple BEC occurs if exactly one eigenvalue of (2.2) is of the order N and
the rest of the eigenvalues are of the order unity.
3. Fragmented BEC occurs if two or more eigenvalues are of the order N and the rest of the eigenvalues are of the order unity.
BEC can be checked depending on the choice of the order parameter defining the condensate. Hence, I include the two different approaches to check the conden-sate.
First criterion for the BEC with an associated order parameter is based on using the boson field operators. The order parameter for this approach can be written as [20];
where ˆψ (~r) is the boson field operator satisfying the commutation relations ∀~r, ~r′; h ˆ ψ (~r) , ˆψ~r′i = hψˆ†(~r) , ˆψ†r~′i= 0 (2.4) h ˆ ψ (~r) , ˆψ†~r′i = δ~r − ~r′ (2.5)
Note that the order parameter, i.e. ˆψ (~r) can be expanded to include time as well. In order for BEC to occur, in the limit ~r − ~r′ → ∞, right hand side of the (2.2)
must have a value different than zero [20].
Second criterion for the BEC with an associated order parameter, on the other hand, based on including density and phase terms. The order parameter for this approach can be written as [20];
ψ (~r) ≡ pN0χ0(~r) (2.6)
Obeying the normalization condition;
Z
d~r|ψ (~r) |2 = N0 (2.7)
where χ0(~r) is a single-particle wavefunction which is orthonormal and can be
defined as [20];
χ0(~r) ≡ |χ0(~r) |eiϕ(~r) (2.8)
To check the BEC in this approach, one can use the single-particle density matrix in the limit |~r − ~r′| → ∞ [20]; lim |~r−~r′|→∞ ρ1 ~r, ~r′= f∗(~r) f~r′+ ˜ρ 1 ~r, ~r′ (2.9) In the limit |~r − ~r′| → ∞, ˜ρ 1
~r, ~r′→ ∞ and f (~r) goes to a non-zero value if
the condensate occurs [20].
2.2
Superfluid Phase
In 1938, two different groups (Kapitsa in Moscow, and Allen Misener in Cam-bridge) simultaneously realized a peculiar behaviour of 4He below the λ-point,
i.e. ≃ 2.17K. They observed that the liquid flowed easily without friction through a narrow channel between two bulk reservoirs [20, 21]. This behavior was later labeled as superfluid by Kapitsa.
The more clear and modern definition of superfluidity, on the other hand, can be understood via observing two conceptually different cases, which can also be called as generalized BEC [20]. These two cases are considered on a multiply connected geometry, e.g. annular region between two concentric cylinders, unlike the narrow channel between two bulk reservoirs.
Due to the fact that superfluid velocity is not a directly measurable quantity, it is suitable to define another quantity in which one can track the existence or the absence of superfluid phase, i.e. current density [20];
~
J (~r) = ρs(~r) ~vs(~r) + ρn(~r) ~vn(~r) (2.10)
where ρn is normal fluid density, whereas ρs superfluid density, and ~vn is normal
This is actually called the two fluid model, and it can be explained with hydro-dynamic equations
Two give a full picture of the phenomena, one needs to define another condition in which phase of the condensate wavefunction is also included. For the spinless case ( it can be written for the spinfull case as well) candensate wavefunction can be written as [20];
χ0(~r, t) ≡ |χ0(~r, t) |eiϕ(~r,t) (2.11)
Hence, the current carried by condensate;
~
J (~r, t) ≡ N0|χ0(~r, t) |2
¯h
m∇ϕ (~r, t)~ (2.12)
where N0 is number of particles in the condensate.
Note that the ratio J(~ρ(~r,t)r,t) has the dimensions of velocity which is defined as
su-perfluid velocity [20];
~vs(~r, t) ≡
¯h
m∇ϕ (~r, t)~ (2.13)
1. Irrotationality: ~ ∇ × ~vs(~r, t) = 0 (2.14) 2. Onsager-Feynman quantization: I Cd~l·~v s(~r, t) = nh m (2.15)
where n is an integer and also called winding number.
Hence, the more clear definition for superfluid phase, accepting the multiply connected geometry, now, can be explained with two cases as mentioned [20];
1. Hess-Fairbank Effect:
While an annulus is being rotated with an angular velocity ω, cool the system down to the λ-point, and wait for the system to reach thermal equi-librium. Since the temperature is above the λ-point, Helium will behave like a normal liquid. Hence, the current density will be ∝ ρ~vn with ~vn= ~ω × ~r,
and the angular momentum is ~L = Iclassical~ω .
Now, cool the system down below the λ-point. One expects to see an-other phase; i.e. He-II [20, 21], the superfluid phase. Due to the Onsager-Feynman quantization condition, (2.15), we have discrete values for vs, i.e.
vs = nωcR, where n is the winding number defined in (2.15).
The value of n in which the system prefers can be determined by minimizing the effective energy [20];
ˆ
Hef f ≡ ˆHlab− ~ω ·L~ˆ (2.16)
where the annulus is stationary in lab frame, which is given by;
ˆ Hlab = X i p2 i 2m + Vext(~ri) + 1 2 X i,j U (|~ri− ~rj|) (2.17)
Minimization of ˆHef f gives an expression depending on the winding num-ber [20]; ε = ρs(~r) R2 1 2n 2ω2− nωω c (2.18) This expression is further minimized with respect to n by choosing n as integer values of the ratio ωωc, i.e. n ≃ intωωc [20]. Contribution of superfluid and normal components to the total angular momentum can be identified depending on which value that ω chooses, i.e. [20];
(a) For ω < ωc
2 ; ~vs = 0 and the total angular momentum is reduced by a
factor ρn
ρ ≡ gn(T ) with L(T ) = gn(T ) Iclassicalω.
(b) For ω > ωc
2 ; ~vs 6= 0 and the total angular momentum is reduced by ρs
ρ ≡ gs(T ) with superfluid contribution to the above equation to the
angular momentum is L(T ) = gs(T ) Iclassicalnωc.
Hence, the total angular momentum is [20];
L(T ) = [gn(T ) ω + gs(T ) nωc] Iclassical (2.19) with n ≡ intω ωc + 1 2 . 2. Metastability of Supercurrents:
This time the annulus is being rotated with a much larger angular velocity ω0 ≡ ω ≫ ωc [20]. Again, cool the system down to the λ-point. vs, by
the Onsager-Feynman constraint, will take quantized values which are close to int ωc
ω
. However, the value of vs, due to ω0 ≫ ωc, should be so small
that any contribution made by the superfluid component will be small and angular momentum would be only L(T ) ≃ Iclassicalω0.
If we further cool the system down below the λ-point, and stop the rotation of the annulus, vn= 0, but superfluid component persist for a while, then,
n ≡ intω0
ωc +
1 2
L(T ) ≃ gs(T ) ω0Iclassical. (2.20)
That is even though it is not the equilibrium one, system has a persistent superfluid circulation which leads to metastable superflow.
In order to determine the superfluid fraction, one needs to consider the system response to the boundary condition. It can, basically, be measured either by cal-culating the free energy change in a periodic system which is based on winding number calculation or by calculating the momentum-density correlation func-tion [5].
1. Winding number Approach:
This approach is based on calculating the density matrix for the moving walls, i.e. rest frame for the walls; ρv [5];
ρv = e−β ˆH, H =ˆ
X
i
(~pi− m~v)2
2m + V (2.21)
where V is the interaction potential.
Response of the fluid to the boundary motion can be written in terms of the total momentum operator, ~P [5];
< ~P >v= ρn ρ N m~v = T rnP ρ~ v o T r {ρv} (2.22) where ρn is the normal component of the fluid.
This equation can also be written with respect to the superfluid fraction [5];
ρs ρ = ∂ (Fv/N ) ∂ 1 2mv2 (2.23)
where Fv is the free energy, and free energy change can be computed as [5]; △Fv N = 1 2mv 2ρs ρ + O v 4. (2.24)
Note that by (2.21), the density matrix satisfies the Bloch equation, hence, obey the periodic boundary condition requirement which brings a phase factor in front of the density matrix coming from the path ending as [5]; eimh¯~v·~L.
This factor introduces the so called winding number, W [5];
N X j=1 ~ r′ j − ~rj = ~W L (2.25)
where ~rj is the initial point and ~r′j is the destination point.
Free energy change can be calculated by using the winding number [5];
e−β△Fv = R d~rρ~v(~r, ~r; β) R d~rρ~v=0(~r, ~r; β) =< eim¯h(~v· ~W L) > (2.26) By using the △Fv
N expansion in (2.24), for small velocities, the above relation
can be written for a d-dimensional system as [5];
ρs ρ = m ¯h2 < W2 > L2−d ρdβ (2.27)
2. Momentum-Density Correlation Function Approach:
This approach is based on the momentum response of the fluid to the bound-ary motion. For a system having periodic boundbound-ary conditions, only the normal component of the fluid responds, and hence, expanding (2.22) to the first order in v gives [5];
ρn
ρ N m~v =< ~P >v= β < ~P ~P > ·~v (2.28) In terms of momentum density, i.e. ~p(~r), (2.28) can be written as [5];
< ~p (~r) >v= β
Z
over all volume
d3r′~v· < ~p (~r) ~pr~′>
v=0 (2.29)
In an isotropic liquid, the normal component of the fluid can be written in terms of momentum density correlation function as [5];
ρn=
β 3m
Z
d3r < ~p (~r) · ~p (0) > (2.30)
3. Single-Particle Delocalization Approach:
Apart from previously used techniques, there is another way to detect the existence of superfluid phase based on Drude weight.
This approach is based on the Drude weight expression which was basically introduced to distinguish metals from insulators, and can be expressed in the form [7, 8]; Dc = π V ∂2E (φ) ∂φ2 φ=0 (2.31) where E is the ground state energy and φ is the phase introduced as per-turbation.
Based on (2.31), a more general relation can be written for the second order derivative of the ground state energy as [10];
∂2E (φ) ∂φ2 φ=0 = i N X j=1 < Ψ|∂kj, ∂xj |Ψ > − lim △X,△K→0 1 △X△K[Im ( ln < Ψ|e i△K ˆXei△X ˆK |Ψ > ei△X ˆK !) + Im ( ln < Ψ|e i△X ˆKe−i△K ˆX|Ψ > ei△X ˆK !) ] (2.33) where ˆK =PNj=1kˆi and ˆX = PN
j=1xˆi, with ˆki and ˆxi are single momentum
and position operators for each particle respectively.
In order to observe the existence of superfluid phase, (2.33) can be used in terms of the sum over expectation values of single momenta, and after taking the limit indicated, one can obtain [9, 10, 22, 23];
σ2x= − 2 (△K)2Re
n
ln< Ψ|e−i△K ˆX
|Ψ >o (2.34) If the wavefunction is an eigenstate of the one-body position shift operator, i.e. e−i△K ˆX, then, this expression contributes to the superfluid weight with
Chapter 3
Bose Hubbard Methods
3.1
Analytical Approach
3.1.1
Mean-Field Theory
Within a correlated system, the motion of each individual particle depends on all the others. To simplify such a system, a physical model is introduced in which correlations between the particles are not entirely included, but, instead, they are included on average. Hence, the effect of the other particles is introduced as mean-field, and the model can be treated as a single particle model.
Mean-field Hamiltonian, i.e. ˆHM F, for the Bose-Hubbard model can be
intro-duced as [15]; ˆ HM F = X i −ψ∗Bˆbi− ψBˆb†i +U 2 X i ˆ ni(ˆni− 1) − µ X i ˆ ni (3.1)
where ψB is a variational operator. ˆb†i , ˆbi are boson creation and annihilation
operators respectively. ˆni = ˆb†iˆbi is the number operator which gives the number
potential.
The variational operators in (3.1) represent the neighboring effects, and they break the U(1) symmetry; hence, phase due to the broken-symmetry is introduced as superfluid phase [15].
In order to approximate the ground state energy of the BHM, a suitable ground state wavefunction must be introduced with an optimum value for the variational operator which minimizes the ground state energy.
Since we can consider the model as a single particle problem, a suitable ground state wavefunction is simply the product of single-site wavefunctions, and hence the ground state energy is [15];
E0 L = EM F(ψB) L − zJ < ˆb †>< ˆb > + < ˆb > ψ∗ B+ < ˆb† > ψB (3.2)
where z is the coordination number, and L is the lattice site, with an optimum value ψB= zJ < ˆb > [15].
Basically, three limits can be applied to the Hamiltonian;
1. J=0 (ψB=0):
For this case, since the variational operator is excluded, sites become decou-pled which gives exact result for the MFT [15]. The Hamiltonian contains only the number operator, thus, the problem is reduced to find only these boson occupation numbers which minimizes the Hamiltonian.
Since ˆn is a good quantum number, ground state wavefunction can be rep-resented by these occupation numbers |mi = n0 Uµ
n0 µ U = 0, Uµ < 0 1, 0 < Uµ < 1 ... n, n − 1 < Uµ < n (3.3) 2. J6= 0 but small:
For this case, since perturbation commutes with the Hamiltonian, system will evolve having exactly the same eigenvalues with the adiabatic increase of J. Hence, the exact result would be [15];
< ˆb†iˆbi >= n0
µ U
(3.4) This result is responsible for the island existing in the phase diagram which are called as Mott insulators [15]. This phase is incompressible with∂<N >
∂µ =
0.
3. J6= 0 (ψB 6=0):
For this case, ground state is delocalized over the lattice. Hence, one cannot use the same ground stated as indicated in (3.3), but, instead, can use [18];
|ΨN >U=0 = 1 √ N ! √1 L! X − → R ˆ a†−→ R N |0 > (3.5)
Since the particle density does not take quantized values, it can change with µ, i.e. ∂<N >
∂µ 6= 0 which defines compressibility [15].
To determine the phase boundaries; as the variational operator, ψ is increased continuously, numerical analysis shows that the Mott insulator phase is a sec-ond order phase transition which can be explained using Landau theory [15]. Expansion of the ground state energy with respect to ψ gives [15];
E0 = E00+ r|ψ|2+ O |ψ|4
∗
(3.6)
where the coefficient r can be found by using the second order perturbation the-ory [15]; r = Γ0(1 − zJΓ0) (3.7) where Γ0 is Γ0 = n0 Uµ + 1 U n0 Uµ − µ + n0 Uµ µ − U n0 Uµ − 1 ! (3.8)
For r=0, phase boundary is found.
3.1.2
Perturbative Methods
The perturbative method is another way to treat the Bose-Hubbard Hamilto-nian. With this method, one can use either strong or weak coupling approaches. In this thesis, strong coupling approach is covered based on the calculations of Freericks [25].
For the strong coupling limit, kinetic energy vanishes and each site has a fixed number of bosons; n0.
For such a system, let the chemical potential to be parametrized as [25]; µ = (n0+ δ) U where n0 is ground state boson occupancy, and δ is the deviation from
integer filling.
determined by calculating the energy of the Mott insulator phase and the de-fect phase, which occurs due to adding hole (δ > 0) or particle (δ < 0) to the system [25], and then, treating the kinetic energy term perturbatively.
Based on the two cases for the defect phase, i.e. δ > 0 and δ < 0, one needs to calculate three different relations: EM I, Edef ectδ>0 , Edef ectδ<0 to determine the
phase boundary.
Relevant wavefunctions which are to the zeroth order in t
U for calculating the
energy relations are given as [25];
1. Wavefunction which belongs to the Mott insulator phase: |ΨM ott(n0) >(0)= N Y i=1 1 √ n0! ˆb† i n0 |0 > (3.9)
2. Wavefunction which belongs to the particle for defect phase:
|ΨDef(n0) >(0)(δ<0)= 1 √ n0+ 1 X i fiˆb†i|ΨM ott(n0) >(0) (3.10)
3. Wavefunction which belongs to the particle for defect insulator phase: |ΨDef(n0) >(0)(δ>0)= 1 √n 0 X i fiˆbi|ΨM ott(n0) >(0) (3.11)
where N is the number of sites in lattice, fi is the eigenstate of the hopping matrix
tij with the lowest eigenvalue.
Energy differences to the third order in Ut between the Mott insulator and defect phases are given as [25];
1. For an extra particle:
EDef(δ<0)(n0) − EM ott(n0) = −δ(particle)U − zt (n0+ 1)
+ zt 2 U n0(5n0+ 4) 2 − z2t2 U n0(n0+ 1) + t 3 U2n0(n0+ 1) [ −2z3+25 4 z 2 − 4z n0 + −z3 +7 2z 2 − 2z ] (3.12) 2. For an extra hole:
EDef(δ>0)(n0) − EM ott(n0) = δ(hole)U − ztn0
+ zt 2 U (n0+ 1)(5n0+ 1) 2 − z2t2 U n0(n0+ 1) + t 3 U2n0(n0+ 1) [ −2z3+ 25 4 z 2 − 4z n0 + −z3+11 4 z 2 − 2z ] (3.13) where EM ott(n0) = N −δUn0− 1 2U n0(n0 + 1) − zt2 U n0(n0+ 1) (3.14)
Phase boundary between the Mott insulator phase and the superfluid phase can be found by setting the energy difference to zero where these two branches meet at [25]; δ(particle)(n
0) + 1 = δ(hole)(n0).
Hence, upper and lower boundaries of the Mott insulator lobe in one dimension are given as [25];
1. Upper boundary: δ(particle) n0, t U = −2 (n0+ 1) t U + n20 t U 2 + n0(n0+ 1) (n0 + 2) t U 3 (3.15) 2. Lower boundary: δ(hole) n0, t U = 2n0 t U − (n0+ 1)2 t U 2 + n0(n0+ 1) (n0− 1) t U 3 (3.16)
A similar approach can be done for weak coupling limit.
3.2
Computational Approach
Dealing with the physics of many-body systems consisting of a large number of interacting particles is in general difficult. Finding an exact solution for such systems as the dimensionality and the total particle number increases becomes impossible. Therefore, a suitable computational approach to such systems is necessary.
In this section, most common computational approaches to the many-body sys-tems are discussed.
3.2.1
Exact Diagonalisation Method
In order to find the eigenvalues of a n-dimensional many-body Hamiltonian one needs to solve a characteristic polynomial with degree n in which finding an exact solution is not possible for n>4 [26].
To be more quantitative, for Hubbard model, in general, with N particles there are 4N states. This brings a limitation on the lattice site as it costs computer
time and memory. For instance for the ultra-cold atom systems L ∼= 22− 25, for square lattice at half-filling it is L ∼= 20, and for the triangular lattice L ∼= 21 [27]. Thus the exponential growth of the matrix ˆH even with small lattice sites makes it hard to calculate the eigenvalues with the usual diagonalisation methods. One suggestion to this problem is that by using the symmetries of the model, one can find a unitary transformation which has the same characteristic polynomial, i.e.
ˆ
H → U†HU , and find its eigenvalues instead.ˆ
It is necessary to construct the U matrix in an iterative way until the matrix ˆH becomes diagonal, i.e. ˆH → U1†HUˆ 1 → U2†U1†HUˆ 1U2 → ... In order to diagonalise
the Hamiltonian Lanczos type algorithms can be used. To do this, one should choose a convenient basis function first and then recursively produce new states until the ground state energy is converged.
As for the low temperature systems, the most relevant eigenstates are either the ground state or the lowest lying excited states. Thus, the initial random choice of state can be chosen with a finite overlap with the ground state. That is [28]
|ψm+1 >= ˆH|ψm > −αm|ψm > −βm2|ψm−1 > (3.17) with coefficients; αm = < ψm| ˆH|ψm > < ψm|ψm> βm2 = < ψm|ψm > < ψm−1|ψm−1 > (3.18)
Then, as a final step, diagonalise the obtained sparse matrix.
With the exact diagonalisation method, one can calculate static quantities like correlation functions or dynamical quantities like density of states.
3.2.2
Quantum Monte Carlo Method
Dealing with quantum systems is more difficult than dealing with the classical systems because in the former case one does not know the exact distribution which is to be sampled however one knows the exact Hamiltonian to be solved. QMC methods are generally based on a random walk process. Simulation starts with a random but reasonable configuration of the system. Then, probability of each configuration is calculated based on Metropolis algorithm, which will be covered in Chapter-4 in detail. With the help of this algorithm, one extracts the ’good’ probabilities, i.e. accepted ones which lead convergent expectation values, and is able to calculate the mean values for the physical system.
Except for the projection method algorithm, which will be covered in Chapter-4 in detail, the most common algorithms which constitute the basis of QMC simu-lations are ’Discrete-time world-line algorithm’ and ’Stochastic series expansion algorithm’.
3.2.2.1 Discrete-Time World-Line Algorithm
This algorithm is based on the path integral formulation of the partition function in imaginary time. The aim is to compute the physical observables either in the canonical ensemble or in the grand canonical ensemble.
The core idea is to use the Suzuki-Trotter decomposition. For the case of a 1 dimensional system, when only nearest neighboring site coupling is allowed, i.e. h
ˆ
Hi,i+1, ˆHj,j+1
i
= 0 where j > i + 1. Thus, one can split ˆH into even and odd terms like ˆH = ˆHeven+ ˆHodd and expand the partition function as [29, 30];
Z ≃ T r ( L Y n=1 e−∆τ ˆHevene−∆τ ˆHodd ) = X n1,...,n2L < n1|exp −∆τ ˆHeven |n2L>< n2L|exp −∆τ ˆHodd |n2L−1 > ... < n3|exp −∆τ ˆHeven |n2 >< n2|exp −∆τ ˆHodd |n1 > (3.19) where L is the lattice site, ∆τ = βL is the imaginary time, and {|ni >} is the
complete basis set for each imaginary time interval.
This decomposition leads to a checkerboard picture of space-time to track the movements of particles along the world-lines which are moved via local up-dates [31]. However, the Suzuki-Trotter decomposition brings an error term on the order O(∆τ2). In order to overcome this error, one needs to introduce the
continuous time limit, i.e. ∆τ → 0, [32, 33];
Z = T rne−β ˆHo = T rne−β ˆHDe− Rβ 0 dτ ˆHOD(τ ) o = T r e−β ˆHD 1 − Z β 0 dτ ˆHOD(τ ) + 1 2 Z β 0 dτ1 Z β τ1 dτ2HˆOD(τ1) ˆHOD(τ2) + ... (3.20) where D stands for diagonal and OD stands for off-diagonal with β = k1
BT.
Note that in the interaction representation, one can write the time-dependent off-diagonal Hamiltonian as ˆHOD(τ ) = eτ ˆHDHˆODe−τ ˆHD.
However, still, such local updates on the checkerboard picture do not change the global properties like number of world lines as a cost of using canonical ensemble. In order to use grand canonical ensemble, global updates must be introduced with the so called loop algorithm [34] and its continuous time limit version [35].
3.2.2.2 Stochastic Series Expansion
This algorithm is based on power-series expansion of the partition function [36];
Z = T rne−β ˆHo = ∞ X n=0 βn n!T r − ˆHn = ∞ X n=0 βn n! X {m1,...,mn} X {b1,...,bn} < m1| − ˆHb1|m2 >< m2| − ˆHb2|m3 > ... < mn| − ˆHbn|m1 > (3.21)
where b is the bond index. Note that one can obtain (3.5) by setting ˆHD = 0
and ˆHOD = ˆH in (3.4).
As it is seen from (3.4) and (3.5) that in world-line algorithm only the off-diagonal is treated as perturbation series whereas in stochastic series expansion whole Hamiltonian is treated as perturbation series. In practice SSE is preferred because in continuous world-line algorithms, one has to deal with the high-precision values of the imaginary time [37].
3.2.3
Density Matrix Renormalization Group Method
DMRG, which was developed by S. R. White [38,39], is based on an iterative and also a variational method in which only the most significant states are considered. The aim is to divide the system into blocks and treat every block separately. While doing this, the interactions among the blocks must also be considered. Starting point for DMRG is the block renormalization which follows the steps basically [40];
1. Create an initial block-A, i.e. ˆHA with length l acting on an n-dimensional
Hilbert space.
2. Then, create a compound block-AA, i.e. ˆHAA with length 2l. The
com-pound block Hamiltonian consists of two block Hamiltonians which has dimensionality = n2.
3. Diagonalise ˆHAA and find the n lowest-lying eigenvectors.
4. Project ˆHAA on the truncated space which is spanned by n lowest-lying
eigenvectors, i.e. ˆHAA → ˆH
′
AA
5. Start from 2l → l and ˆH′
AA → ˆHA. Till the lattice site is reached.
In this method, blocks are considered as independent systems and each block has its own boundary condition. When these blocks are combined same boundary conditions for each block cannot be applied this time as it will not give the true ground state [41]. Thus, instead of just calculating the ground state of the block itself, one should consider the ground state of the compound system and the environment, i.e. block, and after finding the ground state of the super-block, found state is mapped on the block and the block space is truncated with the following formulation [41];
Compound system state can be written as;
|Ψ >=X
i,j
λij|αi > ⊗|βj > (3.22)
where |αi > is a state in block space and |βj > is a state in environment space.
As a next step, environment is traced out and the density matrix for the block only is calculated;
< α′i|ρBlock|αi >=
X
j
Note that the density matrix, i.e. ρBlock, must be;
1. Self-adjoint: ρ†= ρ,
2. Semi-positive definite: ρ ≥ 0, 3. Tracing out to unity: T r {ρ} = 1.
The goal is to find the states which have highest eigenvalues so that the ground state of the super-block is described properly.
Clear definition of steps can be summarized as;
Define m ≡ number of states in a block and n ≡ number of states on a site. 1. Introduce left and right Hamiltonians; i.e. ˆHL and ˆHR acting on an
m-dimensional Hilbert space.
2. Introduce the interactions as left-center and right-center; i.e. ˆHLCand ˆHRC.
3. Then, introduce the super-block Hamiltonian; ˆHSB which is consist of ˆHL,
ˆ
HR, ˆHLC, and ˆHRC with dimensionality = m2n2.
4. Diagonalise ˆHSB and find the ground state.
5. Calculate ρBlock for left and also for right part.
6. Calculate the m-eigenvectors having the highest eigenvalues for left and also for right part.
7. Map left part Hamiltonian; i.e. HˆL, ˆHLC. on m-dimensional truncated
space spanned by the eigenstates which are found in step-6, and do it for the right part of Hamiltonian as well.
There are two types of DMRG; one being the infinite-size DMRG and the other is finite-size DMRG. For further details one can check [41].
DMRG is one of the good and powerful technique for one dimensional problems as the method is based on a low-entanglement approximation. Thus, it is usually preferred to obtain exact solutions in one dimensional.
Chapter 4
Algorithm
4.1
Variational Monte Carlo Method
Quantum Monte Carlo methods allow to calculate expectation values with the help of stochastic sampling by using the so called Metropolis algorithm [42]. This algorithm generates Markov chains, i.e. random walks, over a configuration space. Each configuration is sampled based on a stationary probability distribution. VMC, on the other hand, is one of the QMC methods in which Metropolis al-gorithm is directly used to describe the ground state properties of the system stochastically based on a suitable trial wavefunction.
As a historical side note, VMC method was first applied to a bosonic many-body system to observe the ground state properties of the4He [43], and, it was applied
to the Hubbard model [44] by introducing the celebrated Gutzwiller wavefunction, and the final basic contribution is done by using the square of an anti-symmetric wavefunction to sample the configuration space as an equivalent approach for fermionic many-body problems rather than bosonic ones [45, 46].
4.1.1
Metropolis-Hastings Algorithm
The basic idea of this algorithm is to sample over the configuration space based on acceptance or rejection criteria but keep only the good samples. Applying this idea to the VMC leads to an algorithm based on the following steps;
1. Choose a set of coordinates in the Markov chain; {xi}j randomly (or from
the obtained set from previous configuration).
2. Then, suggest a move with a trial set of coordinates; x′i
j
. The probability of accepting the move is
P = min 1, ψx′ i j ψ {xi}j 2 (4.1) where ψ {xi} j
is the variational wavefunction of the system with a configu-ration {xi}
j
.
3. Generate a random number r st. 0 < r ≤ 1. 4. If ψnx′ioj ψ{xi}j 2
≥ r, accept the move: {xi}j+1 =
x′i
j
. Else, reject the move: {xi}j+1= {xi}j.
5. Then, suggest a new move (step-2) and repeat the process.
Based on the central limit theorem, for large enough samplings, average quantities calculated with the Metropolis algorithm give reliable estimates of the expectation values.
4.2
Adaptation of the VMC Method to the
Bose-Hubbard Model
The following derivations can be applied to any dimension; however, in this thesis results are obtained for a one dimensional system only.
We adapated VMC to the BHM with two different variational wavefunctions, namely; the Baeriswyl and Gutzwiller-Baeriswyl variational wavefunctions.
4.2.1
Baeriswyl Projection
Baeriswyl variational wavefunction is defined as;
|ΨB >= e−α ˆT|ψU=∞> (4.2)
where α is the variational parameter and ˆT = −tP<i,j>ˆb†iˆbj + h.c.
with t as the hopping parameter between nearest neighbouring sites; < i, j > on real space. The operator ˆT can be also represented on momentum space as; Tˆ = −tPLk=02 cos 2πLk
with L being the lattice size.
The term e−α ˆT is called the Baeriswyl projection operator and we projected it on
U = ∞ state, i.e. over a localized state [47].
4.2.2
Gutzwiller-Baeriswyl Projection
The Gutzwiller-Baeriswyl variational wavefunction is defined as;
where γ is the second variational parameter, and ˆN = PLi=0ni(ni− 1) with
ni = ˆb†iˆbi on real space.
The term e−γ ˆN is called the Gutzwiller projection operator for a bosonic
sys-tem [48].
This wavefunction is similar to the one which Otsuka suggested for fermionic systems [49];
|ΨGB >= e−α ˆTe−γ ˆN|φ > (4.4)
where in this case φ is the non-interacting Fermi sea.
Note that e−γ = 0 corresponds to the insulating state for fermions [50].
4.2.3
Calculation of The Observables
In order to calculate the observables, we introduce three different coordinates as left, i.e. xL, center , i.e. xC, right , i.e. xR, for each single particle.
These coordinates form a complete set, i.e.
L
X
i=0
|xL >< xL| = 1 (4.5)
where L is the lattice site, and same condition also applies for center and right coordinates.
As expectation value, in the most general form, can be written as;
< ˆA >= < Ψtr| ˆA|Ψtr > < Ψtr|Ψtr >
where Ψtr is variational wavefunction.
1. The expectation value for the Baeriswyl projection: The operator ˆA is considered as diagonal in real space.
< ˆA > = P L,C,R< Ψ∞|xL >< xL|e−α ˆT|xC >< xC| ˆA|xC >< xc|e−α ˆT|xR >< xR|Ψ∞> P L,C,R< Ψ∞|xL >< xL|e−α ˆT|xC >< xc|e−α ˆT|xR >< xR|Ψ∞ > = P L,C,RP (xL, xC, xR) ˆA (xC) P L,C,RP (xL, xC, xR) (4.7) where L, C, and R stands for left, center, and right coordinates.
The probability function for an accepted move in (4.7) can be written as;
P (xL, xC, xR) = Ψ∞(xL) Ψ∞(xR) K (|xL− xC|) K (|xC− xR|) (4.8)
where Ψ∞(xL) =< Ψ∞|xL>, and the propagator K is defined as;
K|x − x′|=< x|e−α ˆT|x′ >= L Y i=1 < xi|e−α ˆT (1) |x′i > (4.9)
with ˆT(1) being the single particle operator.
By using the fact that the operator ˆT(1) is diagonal in momentum space,
one can write propagator in this space as;
K|x − x′|= L Y i=1 L X k=1 1 Le −αǫkeik x′i−xi (4.10) where ǫkn = −2t cos 2πn L x with n=0,..,L-1.
2. Expectation value for the Gutzwiller-Baeriswyl projection: The operator ˆA is considered as diagonal in real space.
< ˆA > = P L,C,RΨ∞(xL)K (|xL− xC|) e−γ ˆN(xC)A(xˆ C)e−γ ˆ N(xC)K (|x C − xR|) Ψ∞(xR) P L,C,RΨ∞(xL) < xL|e−α ˆT|xC >< xC|e−2γ ˆN|xC >< xC|e−α ˆT|xR > Ψ∞(xR) = P L,C,RP (xe L, xC, xR) ˆA (xC) P L,C,RP (xe L, xC, xR) (4.11) where eP (xL, xC, xR) = P (xL, xC, xR) e−2γ ˆ N(xC).
Note that dealing with expectation values introduces the products of exponential operators which is similar to the Suzuki-Trotter decomposition of e−τ ˆH, where
ˆ
H = ˆT + ˆV with kinetic and potential energies respectively, which is seen in QMC simulations [46].
As a cost of using variational theory, one needs to minimize the energy relation with respect to the variational parameter(s) so that we can obtain a configuration which is close to the exact ground state, i.e. < ˆH0 > ≤ < ˆH0 >tr, where <
ˆ
H0 > stands for the exact ground state and < ˆH0 >tr stands for the approximated
ground state with a suitable variational wavefunction.
4.2.3.1 Calculation of Kinetic Energy, Potential Energy and Super-fluid Density
In this section, calculations of energies and superfluid density are demonstrated only for the Baeriswyl projection. For the Gutzwiller-Baeriswyl projection, sim-ilar calculations can be done on potential energy and superfluid density. Kinetic energy relation with the Gutzwiller-Baeriswyl projection, on the other hand, is different with Gutzwiller factor which will be demonstrated later.
1. Potential Energy Calculation:
On-site term of the Bose-Hubbard Hamiltonian is ˆV = U2 PLi=1nˆi( ˆni− 1).
Note that potential energy is already diagonal in real space; hence,
ˆ V = U 2 P L,C,RP (xL, xC, xR) ˆni(xC) [ˆni(xC) − 1] P L,C,RP (xL, xC, xR) (4.12) 2. Superfluid Density Calculation:
Based on the last idea used in Section-2, we can calculate the superfluid fraction. We need to calculate the square root of the single particle spread function, which is also called single particle delocalization [51]. Spread func-tion is
< ˆσ2 >= −2 (∆K)2Re
n
ln < e−i∆K ˆX >o (4.13)
where ∆K = 2πL is a shift in momentum, and ˆX represents one-body posi-tion shift operator.
Re-express the momentum shift operator for a more computationally ori-ented way with the Euler’s formula gives
< ei∆K ˆX >=< cos∆K ˆX> +i < sin∆K ˆX>=< Ceiϕ> (4.14) where C is the magnitude and φ is the argument of the spread function. Hence, < ˆσ2 >= 1 (∆K)2ln h < cos∆K ˆX>2 + < sin∆K ˆX>2i (4.15) with ˆX =PLi=1xi(xC).
Hence the single particle delocalization that contributes to the superfluid weight is σ L = √ < ˆσ2 > (4.16)
3. Kinetic Energy Calculation For Baeriswyl Projection: We calculated the kinetic energy observable in two ways;
(a) By taking derivative with respect to the variational parameter, α;
< ˆT > = − ∂ ∂ (2α)ln < Ψ∞|e −2α ˆT |Ψ∞ > = −12P 1 L,C,RP (xL, xC, xR) {X L,C,R P (xL, xC, xR) ... ... 1 K (|xL− xC|) ∂ ∂αK (|xL− xC|) + 1 K (|xC − xR|) ∂ ∂αK (|xC− xR|) } (4.18) (b) By implementing the kinetic energy term directly;
Define: X <i,j> ˆb† iˆbi ≡ X C,C′ |x′C >< xC| (4.19)
Then, substitute (4.18) with its hermitian conjugate to (4.7); hence,
< ˆT >= P L,C,C′,RP (xL, xC, xR) K(|xL−xC′|) K(|xL−xC|) + K(|x C′−xR|) K(|xC−xR|) P L,C,RP (xL, xC, xR) (4.20)
4. Kinetic Energy Calculation For Gutzwiller-Baeriswyl Projection:
For this case, we cannot use the first approach showed in (4.2.3.1.3.a); because for different site indices, i.e. i 6= j, the operators ˆT and ˆN do not commute, i.e. hT , ˆˆ Ni 6= 0. Hence, we need to use the kinetic energy formulation explained in (4.2.3.1.3.b) by adding the Gutzwiller correction. Substitute (4.18) into (4.11), and note that the operator ˆN is diagonal in real space; < ˆT >= P L,C,C′,RP (xe L, xC, xR) K(|xL−xC′|) K(|xL−xC|) + K(|x C′−xR|) K(|xC−xR|) e−γ ˆN(xC′) e−γ ˆN(xC) P L,C,RP (xe L, xC, xR) (4.21)
Chapter 5
Data Analysis and Conclusion
In this chapter, results of the VMC techniques applied to the one dimen-sional Bose-Hubbard model with projected wavefunctions; namely, Baeriswyl and
Gutzwiller-Baeriswyl, are analysed.
As a requirement of the VMC approach, we need to propose a trial wavefunc-tion in order to analyse the ground state properties of the Bose-Hubbard model. Depending on the variational parameter attained to each trial wavefunction, one obtains a set of ground state energies. Among them, the lowest energy state must be chosen, which is described by the best trial wavefunction. In order to choose such a wavefunction, one needs to optimise the energy as a function of the pa-rameter(s); which can be done via methods; steepest descent, parallel tempering
Monte Carlo, energy variance minimization, or conjugate gradient.
In order to observe the ground state properties of the system with the projected wavefunctions, we, first, minimise the energy considering the different choices of the hopping parameter, t>0. Then, we introduce a hypothetical chemical potential, since we are working on the canonical ensemble, as [4];
Hence, we can and will study the ground state phase diagram of the Bose-Hubbard model.
How we control the number of particles is implemented in the algorithm we developed, which is imposed on the ψ∞, regardless of the hypothetical chemical
potential that we introduced.
To study the ground state phase diagram of the model, we introduced two differ-ent variational wavefunctions. First, we analyse the system with the Baeriswyl projected wavefunction, then, we re-analyse the system with a Gutzwiller correc-tion imposed on the Baeriswyl projected wavefunccorrec-tion.
5.1
Baeriswyl Projection Results
Baeriswyl wavefunction, stated in (4.2), is projected onto a localized state and is supposed to introduce hoppings between sites. These hoppings are controlled with the parameter t, as well as α and U. We analysed a system having lattice sites as L = 20 with particle number N = [1, 60] and then L = 40 with particle number N = [1, 160]. Larger lattice site is considered only for the single particle delocalization function, i.e. (4.15).
In order to observe how single particle level is delocalizing with the Baeriswyl projection, one can check the propagator, i.e. (4.10), behavior as;
0 5 10 15 20 25 30 35 40 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 BWH, L=40, N=20 sites prp alpha=0.01 alpha=0.1 alpha=1 alpha=10 alpha=100
Figure 5.1: Propagator vs lattice site for Baeriswyl projection.
Observe that as α gets large, single particle propagator starts to delocalize. Before presenting the phase diagram for the Baeriswyl projection, we checked indirectly that whether a quantum phase transition occurs before introducing the hypothetical chemical potential, i.e. µ in (5.1). Below, one can see how the total ground state energy of the Baeriswyl wavefunction changes with respect to α when we scan over different choices of the hopping parameter for a system with L = 40;
0 0.5 1 1.5 2 2.5 3 −20 −10 0 10 α E TOT Baeriswyl, N=20, L=40 0 0.5 1 1.5 2 2.5 3 −20 0 20 40 Baeriswyl, N=40, L=40 α E TOT t=1 t=0.05 t=0.15 t=0.25 t=1 t=0.05 t=0.15 t=0.25
Figure 5.2: ET OT vs α for Baeriswyl projection.
In Fig.5.2, even though α is just a variational parameter, it is an indication of the order parameter because each α value defines a specific wavefunction, i.e. (4.2) and (4.3), and thus, we would have an order parameter associated to that specific α value. Based on this logic, we do not see any phase transition neither first nor second order. We observe that a global minimum occurs and migrates for larger values of t.
The phase diagram for the Baeriswyl projection with L = 40 lattice sites can be seen in Fig.5.3.
−10 0 1 2 3 4 5 6 7 0.5 1 1.5 2 2.5 3 3.5
Particle Density vs Chemical Potential, BW
µ ρ =N/L t=0.15 0.050 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 1 2 3 4 5 6 t/U µ /U
Phase Diagram For Baeriswyl, L=40
MOTT ρ=3 MOTT ρ=2 MOTT ρ=1 t=0.15 SUPERFLUID
Figure 5.3: Phase diagram for Baeriswyl projection with L = 40.
In Fig.5.3, on the left hand side, we see that for a specific value of the hopping parameter which is chosen arbitrarily how the phase diagram on the right hand side occurs. Observe that for integer values of the boson density, we have a finite range in µ which does not affect the boson density. This case corresponds to the Mott insulating phase. Whereas, for non-integer values of the boson density, we do not have a constant ratio of ρ which corresponds to the superfluid phase. To see how system size affects the phase diagram, observe Fig.5.4 for a system having L = 20 lattice sites;
0.050 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 1 2 3 4 5 6 t/U µ /U
Phase Diagram For Baeriswyl, L=20
t=0.15 MOTT ρ =3 MOTT ρ =2 MOTT ρ =1 SUPERFLUID
Figure 5.4: Phase diagram for Baeriswyl projection with L = 20.
See that the tip of the Mott insulating lobes are different depending on the system size, which gives better results [19] for the larger system size.
5.2
Gutzwiller-Baeriswyl Projection Results
Gutzwiller-Baeriswyl wavefunction, stated in (4.3), is projected onto a local-ized state and is supposed to introduce hoppings between sites as in the case of Baeriswyl projection. For the Gutzwiller-Baeriswyl projection, we analysed a system having lattice sites as L = 40 with particle number N = [1, 160]. Larger lattice site is considered only for the single particle delocalization function, i.e. (4.15).
Again, before presenting the phase diagram for the Gutzwiller-Baeriswyl projec-tion, we checked indirectly that whether a quantum phase transition occurs before
introducing the hypothetical chemical potential, i.e. µ in (5.1). Below, one can see how the total ground state energy of the Gutzwiller-Baeriswyl wavefunction changes with respect to α when we scan over different choices of the hopping parameter for a system with L = 40;
0 0.5 1 1.5 2 2.5 3 −20 −10 0 10 α E TOT Gutzwiller−Baeriswyl, N=20, L=40 0 0.5 1 1.5 2 2.5 3 −20 0 20 40 α E TOT Gutzwiller−Baeriswyl, N=40, L=40 t=1 t=0.05 t=0.15 t=0.25 t=1 t=0.05 t=0.15 t=0.25
Figure 5.5: ET OT vs α for Gutzwiller-Baeriswyl projection with γ = 0.01.
In Fig.5.5, we again observe that a global minimum occurs and migrates for larger values of α. Hence, we do not have any phase transition, neither first nor second order as in the Baeriswyl results.
The phase diagram for the Gutzwiller-Baeriswyl projection with L = 40 lattice sites can be seen in Fig.5.7.
−10 0 1 2 3 4 5 6 0.5 1 1.5 2 2.5 3 3.5
Particle Density vs Chemical Potential
µ ρ =N/L t=0.15 0.050 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 1 2 3 4 5 6 t/U µ /U
Phase Diagram For Gutzwiller−Baeriswyl, L=40
MOTT ρ=1 MOTT ρ=2 MOTT ρ=3 SUPERFLUID t=0.15
Figure 5.6: Phase diagram for Gutzwiller-Baeriswyl projection with L = 40.
In Fig.5.7, as in Fig.5.3, on the left hand side, we see that for a specific value of the hopping parameter which is chosen arbitrarily how the phase diagram on the right hand side occurs. Observe that for integer values of the boson density, we have a finite range in µ which again does not affect the boson density. This case corresponds to the Mott insulating phase. Whereas, for non-integer values of the boson density, we do not have a constant ratio of ρ which corresponds to the superfluid phase.
As a comparison between the two projectors that we applied to the Bose-Hubbard model with same lattice sites, i.e. L = 40, see Fig.5.8;
Figure 5.7: Phase diagram comparison between Gutzwiller and Gutzwiller-Baeriswyl projections with L = 40
The red one with dots belongs to the Baeriswyl results and the blue one with solid lines belongs to the Gutzwiller-Baeriswyl results.
In the PRB paper of Freericks and Monien in 1995, they found the critical value for the Mott lobe as t
U = 0.215 ± 0.01 [19] by applying the QMC techniques.
Even though our results are far beyond what they found, the tips for the Mott lobes are closer to their value for the Gutzwiller-Baeriswyl projection.
5.3
Conclusion
In this chapter, we analysed the results that we obtained from both the Baeriswyl and Gutzwiller-Baeriswyl projections for the one dimensional Bose-Hubbard model even though our method is applicable to any dimension.
Figures 5.2 and 5.5, on the other hand, indicate that for a fixed number of par-ticles, both in commensurate and incommensurate fillings, we do not have phase transition. Hence, we need to introduce another on-site potential, i.e. µ, to control the particle fluctuation and encourage phase transition.
The phase diagrams obtained for both projections give qualitatively good results compared to the results of Freericks and Monien [19] and Scalettar et al. [4]; how-ever, our results are quantitatively weak as the critical t value for the Baeriswyl result is ∼ 0.45; whereas for the Gutzwiller-Baeriswyl it is ∼ 0.5. Although we expect that the Gutzwiller correction to the Baeriswyl projection would improve the critical t value for the tip of the Mott lobes, it gave no good contribution. In order to fix this problem Monte Carlo steps might be increased or different optimization methods, other than minimization of the energy as indicated at the beginning of this chapter, can be done on the projected wavefunctions.
Regarding to the quantitative calculations, our results show that variational Monte Carlo approach with Baeriswyl and Gutzwiller-Baeriswyl projections do give rough results rather than exact. However, the chosen lattice size, which is at most L = 40 for general expectation values might affect the results. As we can observe from the phase diagrams of different lattice sizes, i.e. Fig.5.3 and Fig.5.4, critical value of t changes, and as the lattice size gets bigger, it gives more accurate results regarding to the result in [19].
Future work might include trying these projected wavefunctions with other Monte Carlo methods; like diffusion Monte Carlo or path integral Monte Carlo, which is more common, and see if it gives better results. Preferentially, if variational Monte Carlo is going to be used as a method, one can check how dimensionality, larger Monte Carlo steps, or larger system size affects the results compared to the ones we obtained in this thesis.
Bibliography
[1] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, “Cold bosonic atoms in optical lattices,” Phys. Rev. Lett., vol. 81, pp. 3108–3111, Oct 1998.
[2] A. van Oudenaarden and J. E. Mooij, “One-dimensional mott insulator formed by quantum vortices in josephson junction arrays,” Phys. Rev. Lett., vol. 76, pp. 4947–4950, Jun 1996.
[3] R. Baltin and K.-H. Wagenblast, “Quantum phase transitions for bosons in one dimension,” Europhys. Lett., vol. 39, no. 1, pp. 7–12, 1997.
[4] R. T. Scalettar, G. Batrouni, P. Denteneer, F. H c bert, A. Muramatsu, M. Rigol, V. Rousseau, and M. Troyer, “Equilibrium and dynamical prop-erties of the boson hubbard model in one dimension,” Journal of Low
Tem-perature Physics, vol. 140, no. 3-4, pp. 313–332, 2005.
[5] E. L. Pollock and D. M. Ceperley, “Path-integral computation of superfluid densities,” Phys. Rev. B, vol. 36, pp. 8343–8352, Dec 1987.
[6] J. Cao, “Winding-number effect in path-integral simulations,” Phys. Rev. E, vol. 49, pp. 882–889, Jan 1994.
[7] W. Kohn, “Theory of the insulating state,” Phys. Rev., vol. 133, pp. A171– A181, Jan 1964.
[8] B. S. Shastry and B. Sutherland, “Twisted boundary conditions and effec-tive mass in heisenberg-ising and hubbard rings,” Phys. Rev. Lett., vol. 65, pp. 243–246, Jul 1990.
[9] B. Hetenyi, “Drude weight, meissner weight, rotational inertia of bosonic superfluids: How are they distinguished?,” Journal of the Physical Society
of Japan, vol. 83, no. 3, p. 034711, 2014.
[10] B. Het´enyi, “dc conductivity as a geometric phase,” Phys. Rev. B, vol. 87, p. 235123, Jun 2013.
[11] J. Hubbard, “Electron correlations in narrow energy bands,” Proceedings of
the Royal Society of London. Series A, Mathematical and Physical Sciences,
vol. 276, no. 1365, pp. pp. 238–257, 1963.
[12] H. A. Gersch and G. C. Knollman, “Quantum cell model for bosons,” Phys.
Rev., vol. 129, pp. 959–967, Jan 1963.
[13] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, “Boson localization and the superfluid-insulator transition,” Phys. Rev. B, vol. 40, pp. 546–570, Jul 1989.
[14] p. a. y. Wen Xiao Gang, title=Quantum Field Theory of Many-Body Sys-tems: From the Origin of Sound to an Origin of Light and Electrons. [15] S. Sachdev, Quantum phase transitions. Cambridge: Cambridge University
Press, second ed. ed., 2011.
[16] J. M. Zhang and R. X. Dong, “Exact diagonalization: the bose–hubbard model as an example,” Europhys. Lett., vol. 31, no. 3, pp. 591–602, 2010. [17] T. Comparin, Numerical study of the trapped and extended Bose-Hubbard
models. PhD thesis, Utrecht University, 2013.
[18] I. Bloch, J. Dalibard, and W. Zwerger, “Many-body physics with ultracold gases,” Rev. Mod. Phys., vol. 80, pp. 885–964, Jul 2008.
[19] J. K. Freericks and H. Monien, “Strong-coupling expansions for the pure and disordered bose-hubbard model,” Phys. Rev. B, vol. 53, pp. 2691–2700, Feb 1996.
[20] A. J. Leggett, “Bose-einstein condensation in the alkali gases: Some funda-mental concepts,” Rev. Mod. Phys., vol. 73, pp. 307–356, Apr 2001.
[21] P. Kapitza, “Viscosity of liquid helium below the -point,” Nature, vol. 141, Jan 1938.
[22] R. Resta, “Quantum-mechanical position operator in extended systems,”
Phys. Rev. Lett., vol. 80, pp. 1800–1803, Mar 1998.
[23] R. Resta and S. Sorella, “Electron localization in the insulating state,” Phys.
Rev. Lett., vol. 82, pp. 370–373, Jan 1999.
[24] D. J. Scalapino, S. R. White, and S. Zhang, “Insulator, metal, or supercon-ductor: The criteria,” Phys. Rev. B, vol. 47, pp. 7995–8007, Apr 1993. [25] H. M. J. K. Freericks, “Phase diagram of the bose hubbard model,” Europhys.
Lett., vol. 26, no. 7, p. 545, 1994.
[26] H. Weimer, “Chapter 1: Exact diagonalization.” This is an electronic docu-ment. http://www.itp.uni-hannover.de/ weimer/teaching/QSim2013/exact-diagonalization.pdf. Date retrieved: May 19, 2014.
[27] K. S. Pati, “Density matrix renormalization group
method.” This is an electronic document.
http://www.jncasr.ac.in/ccms/nqm2007/lecturenotes/6day3nov/dmrg Swapan.ppt. Date retrieved: May 19, 2014. Date last modified: May 20, 2014.
[28] E. Dagotto, “Correlated electrons in high-temperature superconductors,”
Rev. Mod. Phys., vol. 66, pp. 763–840, Jul 1994.
[29] M. Suzuki, “Generalized trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems,” Communications in Mathematical Physics, vol. 51, no. 2, pp. 183– 190, 1976.
[30] H. F. Trotter, “On the product of semi-groups of operators,” Proc. Am.
Math. Soc, vol. 10, no. 4, pp. 545–551, 1959.
[31] J. E. Hirsch, R. L. Sugar, D. J. Scalapino, and R. Blankenbecler, “Monte carlo simulations of one-dimensional fermion systems,” Phys. Rev. B, vol. 26, pp. 5033–5055, Nov 1982.
