《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

平成二十六年度 修士論文

Flux quench in the S = 1/2 XXZ chain(S = 1/2 XXZ 鎖模型における磁束クエンチ)

東京大学大学院理学系研究科物理学専攻
中川裕也(指導教員: 押川正毅)

March 10, 2015

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Abstract
In recent years, nonequilibrium quantum systems have become a major subject of study
in condensed matter. Quantum quenches, abrupt changes of Hamiltonian of a system,
oﬀer simple protocols to tackle problems in nonequilibrium situations.
In this thesis, we study a flux quench problem in the spin-1/2 XXZ chain. The flux
quench is a quantum quench where the flux ϕ piercing the XXZ chain is turned oﬀ at
t = 0 suddenly. If we formulate the XXZ chain as a spinless fermion model, the flux ϕcorresponds to a vector potential on each bond and this flux quench can be viewed as
imposing a pulse (delta function) of electric field. Therefore the flux quench generates
some particle (or spin) current at time t = 0. We focus on the time-evolution of the spincurrent after the quench, and calculate it numerically by the infinite time-evolving block
decimation (iTEBD) method.
We find that the long-time limit (t → ∞) and the dynamics of the spin current
depend strongly on the anisotropy parameter ∆ of the XXZ chain and the amount of flux
initially inserted. The long-time limit of the current matches with predictions of linear
response theory as the initial flux decreases. However, the deviation from linear response
theory is largely aﬀected by the sign of the anisotropy of the XXZ chain. Regarding
the dynamics, the current decays in time after the quench, and its relaxation time isdependent on the anisotropy and the initial flux. Furthermore, for large initial flux andanti-ferromagnetic anisotropy, the current oscillates in long time scale. We numericallyfind that the frequency of the oscillation is proportional to |∆|. Remarkably, the dynamics
of momentum distribution of spinless fermions reveals that this oscillation of current isgoverned by excitations deep inside the Fermi sea. This mechanism of oscillations cannot
be captured by the eﬀective Luttinger model corresponding to the microscopic XXZ chain,
which is in contrast with the previous studies on diﬀerent types of quench in the samemodel as ours.

1《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Contents
1 Introduction
1.1 Nonequilibrium physics . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Quantum quench . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 One-dimensional quantum systems . . . . . . . . . . . . . . . . . . . . . .
1.3 Purpose of this study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
5
7
8
8

2 Definition of the problem: Flux quench in the S = 1/2 XXZ chain
2.1 S = 1/2 XXZ chain model . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 XXZ chain with flux . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Flux quench in the S = 1/2 XXZ chain . . . . . . . . . . . . . . . . . . . .
2.2.1 Analysis by linear response theory and the Drude weight . . . . . .
2.3 Previous studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Previous studies on the long-time limit of the current . . . . . . . .
2.3.2 Bosonic flux quench . . . . . . . . . . . . . . . . . . . . . . . . . . .

9
9
11
13
15
15
16
17

3 Numerical Methods
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Tensor networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Basics of tensor networks . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Entanglement entropy and area law . . . . . . . . . . . . . . . . . .
3.2.3 Ground states are well expressed by tensor networks . . . . . . . .
3.3 Matrix product states (MPS) . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Schmidt decomposition . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Construction of MPS for arbitrary states in a system . . . . . . . .
3.3.3 Canonical condition for MPS . . . . . . . . . . . . . . . . . . . . .
3.3.4 Calculation of physical quantity using MPS . . . . . . . . . . . . .
3.3.5 Example of MPS . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Time-Evolving Block Decimation method . . . . . . . . . . . . . . . . . . .
3.4.1 Suzuki-Trotter decomposition for time-evolution operator . . . . . .
3.4.2 Updating canonical MPS after the action of local operator . . . . .
3.4.3 TEBD algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 iTEBD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.1 infinite MPS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.2 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18
18
19
19
21
22
23
24
26
29
31
32
34
34
35
36
38
38
39

2

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

4 Results and Discussions
4.1 Numerical data of the dynamics . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Long-time limit of the current . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Frequency of oscillation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Equation of motion of spin current . . . . . . . . . . . . . . . . . .
4.3.2 Dynamics of momentum distribution nq . . . . . . . . . . . . . . .
4.4 Relaxation time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41
41
45
46
47
47
50

5 Conclusion, Future Work52

A Technical details on numerical calculations
53
A.1 Details on numerical calculations, error estimate . . . . . . . . . . . . . . . 53
A.2 iTEBD algorithm exploiting a symmetry . . . . . . . . . . . . . . . . . . . 55

3

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 1
Introduction
It is of great importance to study nonequilibrium physics, from the viewpoint of puretheoretical interest as well as practical realizations of quantum devices such as quantum
computers. In this chapter, we review background of our study. Nonequilibrium physics
is overviewed and quantum quench problems are introduced in the first section. In thesecond section, we sketch the physics in one dimension. Emphasis is put on a drastic
eﬀect of interactions in one-dimensional quantum systems and a wide variety of methods
to study them. An aim of our study and outline of this thesis are presented in the lastpart of this chapter.

1.1

Nonequilibrium physics

Equilibrium states are defined as states that have become stationary after long timeevolution in closed systems. They are designated by only a few parameters, such astemperature, pressure and volume. In the language of quantum statistical mechanics,
equilibrium states are mixed states (Gibbs states) written asρ̂βGibbs =e−β Ĥ
,Tr(e−β Ĥ )(1.1)

where Ĥ is Hamiltonian of a system and β is inverse temperature. In order to reveal
properties of the equilibrium states, many theoretical methods are proposed: mean-fieldapproximation, renormalization group, Landau-Ginzburg theory, and so on. These methods can capture essential points of experimental data on complicated real materials. So
far, we have much established understanding of the equilibrium states.
On the other hand, less are known for nonequilibrium states. It is partly because
most theoretical methods for equilibrium states cannot be applied to nonequilibrium
situations directly and we do not have enough tools to tackle them especially beyond
linear response regime. Recently, substantial experimental developments in nonequilibrium quantum systems have allowed one to compare theories with experiments at highaccuracy and controllability. For example, optical lattices in cold atom systems [1] cansimulate theoretical models, such as Hubbard model, with broad tunability of parametersand good time-resolutions of the dynamics. As a consequence, nonequilibrium quantum
systems have become a major subject of study in condensed matter physics [2].
4

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 1. Introduction
The field of nonequilibrium physics is really vast. There are various ways to drive
systems into nonequilibrium: coupling them to thermal baths, adding interactions withnoisy environments, applying external fields, etc. Among these protocols, in this study
we choose one called quantum quench.

1.1.1

Quantum quench

Quantum quench, a hot topic in recent studies, is a sudden change of parameter(s) ina quantum system. It is one of the simplest protocols to make systems nonequilibrium
and study the nature of such situations. In typical setup, initial state is prepared as theground state of pre-quench Hamiltonian, and then a parameter of the system is changedin time. The change of parameter induces non-trivial dynamics of the system, and one
studies the dynamics itself and/or its long-time limit.
For example, consider quantum quench of magnetic field in a quantum spin system.
Initial state is taken to be the ground state of Hamiltonian Ĥ(h) in the presence ofmagnetic field h. We assume that the initial state has magnetic order m. Then themagnetic field is quenched to zero at time t = 0 abruptly. Because the ground state ofpre-quench Hamiltonian Ĥ(h) is not an eigenstate of post-quench Hamiltonian Ĥ(h = 0),
non-trivial time-evolution1 is evoked:
Ĥ(h) |ψ(h)⟩ = E0 (h) |ψ(h)⟩ , e−iĤ(h)t |ψ(h)⟩ = e−iE0 (h)t |ψ(h)⟩ ,(1.2)

Ĥ(0) |ψ(h)⟩ ∝/ |ψ(h)⟩ , e−iĤ(0)t |ψ(h)⟩ ̸= e−iE0 (h)t |ψ(h)⟩ .

(1.3)

When the ground state and any thermal state of post-quench Hamiltonian Ĥ(h = 0) do
not exhibit magnetic order, the order parameter m(t) of the state |ψ(t)⟩ = e−iĤ(0)t |ψ(h)⟩is expected to decay to zero if we assume thermalization. In studies of quantum quench,
this kind of equilibration and resulting thermalization are investigated as well as the
dynamics after the quench.
By quantum quench, we can study various kinds of phenomena. One example isuniversal dynamical scaling across a quantum critical point [3]. Universal scaling laws arefound to hold when systems are quenched across quantum critical points, as in the KibbleZurek mechanism for classical systems. Some quantities (density of defects, adiabatic
fidelity, etc.) show scaling laws whose critical exponents are universal and independentof microscopic details of a system. Another example of phenomena studied by quantum
quench is equilibration and thermalization in closed quantum systems. In addition, thedynamics towards (thermal) equilibrium is investigated by employing quantum quench.
Below we review these two in detail.
Equilibration and thermalization in closed quantum systems
Equilibration in closed quantum systems is a phenomenon that quantum states become
stationary in the long-time limit, t → ∞. If such equilibrated states in the long-time limit
look thermal, we call that phenomenon as thermalization. Specifically, thermalization is1

In this thesis, we set ℏ = 1.

5

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 1. Introductiona phenomenon that equilibrated states have the same expectation values for any physicalobservable as those of certain Gibbs states ρβGibbs :∃ βρGibbs ,()
lim ⟨ψ(t)|Â|ψ(t)⟩ = Tr ρβGibbs Â

t→∞for all Â.

(1.4)

We note that the detail of initial state |ψ(0)⟩ is lost in the right hand side and theequilibrated state is specified by only one parameter, inverse temperature β. The initialstate forgets the information about itself during time-evolution and becomes featureless.
Thermalization is rigorously confirmed to happen in some special cases and expected
to occur in generic systems [2,4]. However, for the so-called integrable systems, it is widelybelieved that there is no thermalization. Integrable systems are quantum systems thatare exactly solved by Bethe ansatz, and there exist infinitely many conserved quantities inthe system. Naively, the reason for the absence of thermalization in integrable systems isthat the infinitely many conserved quantities prevent states from forgetting information
about themselves. Nevertheless, for integrable systems we can construct another statistical ensemble to describe equilibrated states after quantum quench: Generalized Gibbs
ensemble (GGE) [5]. GGE is an ensemble that contains additional conserved quantities
other than Hamiltonian,β,{λ}ρ̂GGE =1
ZGGE(
exp −β Ĥ −∑)λi Q̂i

,

(1.5)

i=1

where Q̂i are conserved quantities in a system, λi are Lagrange multipliers and ZGGE isa normalization constant. The point of GGE is that while the total number of states(the dimension of Hilbert space) is exponentially large in system size L, the number ofparameters to designate GGE, β and {λ}, is linearly large in L. In this sense, informationof initial states is lost and thermalization-like behavior appears even in integrable systems.
For some integrable models essentially equivalent to free particles (e.g. transverse
Ising model), it was shown [5] that GGE can represent equilibrated states after quantum
quench. GGE was expected to be valid also in general integrable systems [2]. However,quite recently [6] it was proved that GGE having only local conserved quantities cannot
represent states after quench in a certain model (the S = 1/2 XXZ chain, the same model
as our study). We will review this study (Ref. [6]) in section 2.3 as the previous studies
related to our problem. A general proof (or denial) for GGE is still outstanding, andmany problems remain to be solved in the field.
Dynamics after quantum quench
Equilibration and thermalization phenomena say about only the long-time limit of states.
Transient dynamics towards stationary states after quantum quench is also investigatedby many authors [7,8], although the dynamical behaviors after quantum quench are quite
diverse and it is not easy to conclude universal features of them. In section 2.1, we will
review examples of studies that focused on the dynamics after quantum quench in theS = 1/2 XXZ chain (the model we consider in this thesis).
6

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 1. Introduction

1.2

One-dimensional quantum systems

In this thesis, we study a nonequilibrium problem in one dimension. Here we review thephysics of one-dimensional quantum systems.
One-dimensional quantum systems have been studied for more than 50 years [9]. They
have attracted theorists because of their remarkable features: a strong eﬀect of interactions
on systems and wide availability of theoretical methods and exact results, compared withhigher-dimensional systems. One can regard them as a prototype of strongly-correlated
systems; eﬀects of interactions are quite strong but we can obtain exact (or very detailed)results. Furthermore, they are not only a subject of pure-theoretical studies but also
realistic experimental systems: organic conductors, carbon nanotubes, quantum wires,
optical lattices in cold atomic systems, and so on.
Drastic eﬀect of interactions in one dimension and Tomonaga-Luttinger liquidInteractions play a crucial role in one dimension. The breakdown of Landau’s Fermi
liquid theory manifests the importance of interactions in one dimension. Fermi liquid isa universality class into which a wide variety of metallic systems of interacting fermionsin higher dimensions (> 1) flows under the renormalization group. Fermi liquid theory(as a textbook, we name Ref. [10]) is well established in dimensions higher than one. Its
assumption is the existence of quasiparticle excitations that are adiabatically connected
to excitations in a free system as interactions are turned oﬀ. In one dimension, however,interactions change properties of a system drastically and elementary excitations become
collective and bosonic; they are not connected adiabatically to excitations in a free system.
Consequently, Fermi liquid theory breakdowns in one dimension.
Metallic systems of interacting fermions in one dimension flow into another universality class under the renormalization group: Tomonaga-Luttinger Liquid (TLL). TLL isgapless and critical, and it has no symmetry-breaking orders. It is designated by only two
universal parameters, the renormalized Fermi velocity vF and the Luttinger parameter K.
Low-energy physics of a system belonging to TLL universality class is fully determinedby these two parameters. For example, the exponent of correlation function of fermionsat long distance (low-energy) is written by K: G(r ≫ 1) ∼ 1/r1/2K+K/2 . Furthermore,
one of the most important aspects of TLL is that not only fermionic but also bosonicand spin systems in one dimension mostly exhibit TLL behaviors. TLL is ubiquitous in
one-dimensional quantum physics.
Theoretical methods to study one-dimensional quantum systems
In spite of such strong (non-perturbative) eﬀects of interactions, many analytical andnumerical methods are available to study one-dimensional quantum systems compared tohigher-dimensional systems. We briefly review some of the methods; Bethe ansatz [11]gives exact solutions for one-dimensional Heisenberg model and Hubbard model. The
solvability of these models becomes a branch of mathematical physics under the nameof integrable systems. Bosonization [9] is one of the most powerful theoretical methods
to analyze one-dimensional quantum systems. Peculiarity of one dimension that particles cannot overtake each other allows us to map fermionic systems into bosonic systems
7

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 1. Introduction
and vice versa, which results in making the method applicable to wide classes of models. Conformal field theory [12] predicts behaviors of a system at and nearby a criticalpoint. For one-dimensional systems (1+1 dimension in space-time) conformal invarianceemerging from criticality strongly dictates various physical quantities and we can obtain
exact formulae for such quantities. As for numerical methods, we can employ exact diagonalization (ED) and quantum Monte Carlo for one-dimensional systems [13] in the samemanner for higher-dimensional systems. Moreover, as peculiar methods in one dimension,
density matrix renormalization group (DMRG) [14] and time-evolving block decimation(TEBD) [15] enable us to deal with very large systems (several hundreds of sites) thatcannot be handled by ED or quantum Monte Carlo. We employ infinite-TEBD (iTEBD)
method in this study and it is reviewed in chapter 3.

1.3

Purpose of this study

When we combine two concepts introduced in this chapter, we conceive of a stimulating
field: nonequilibrium physics in one-dimensional quantum systems. Although nonequilibrium phenomena in strongly-correlated systems are formidably diﬃcult to tackle, thewide availability of analytical and numerical methods in one dimension allows us to study
them at depth. The study of nonequilibrium physics in one dimension can be regarded
as a prototype to understand nonequilibrium phenomena in strongly-correlated systems.
As one of the simplest problems in nonequilibrium and strongly-correlated systems, inthis thesis we study a quantum quench problem named flux quench in the S = 1/2 XXZ
chain by exploiting the iTEBD method (numerical method peculiar in one dimension).
The flux quench is a quantum quench of flux piercing a ring of the XXZ chain and it
induces spin current in the system. We study the long-time limit and dynamics of theinduced spin current after the quench.

1.4

Outline of the thesis

This theses is organized as follows. In chapter 2, the S = 1/2 XXZ chain model is reviewed
and the flux quench problem is defined. Besides, we introduce the previous studies related
to our problem. In chapter 3, numerical methods employed in this study are reviewed. It
is almost independent of other contents in this thesis, so those who are not interested indetails of numerical calculations can skip it. It also can be read as an independent reviewof the TEBD method. In chapter 4, the results of numerical calculations are shown.
Discussions on the numerical data are also presented. In chapter 5, we state a conclusionof this study and future works are discussed. Technical details of numerical calculations
are described in appendix A.

8《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 2
Definition of the problem: Flux
quench in the S = 1/2 XXZ chain
In this chapter, we review properties of the S = 1/2 XXZ chain model and define a
problem of our study, flux quench.

2.1

S = 1/2 XXZ chain model

The S = 1/2 XXZ chain model is a generalization of Heisenberg model which describes
quantum magnets in one dimension. Hamiltonian of this model is defined as
HXXZ = −
where

N
−1
∑

N
−1
∑
( x x
)
y y
z z
Si Si+1 + Si Si+1 + ∆Si Si+1 = −

i=0

i=0

(

)
)
1( + −− +z z
S S + Si Si+1 + ∆Si Si+1
2 i i+1
(2.1)

)
)
)
(
(
(
1 0 −i
1 1 0
0 1
y
z
, Si =, Si =, Si± = Six ± iSiy
1
0
i
0
0
−1
2
2
2

1
Six =

(2.2)

x,y,z
=
are spin-1/2 operators acting on site i. We impose periodic boundary condition SN
S0x,y,z . ∆ is called the anisotropy of the chain. ∆ > 0 (∆ < 0) means ferromagnetic(anti-ferromagnetic) interactions between neighboring sites. When ∆ is ±∞, the model
is equivalent to the Ising model.
This Hamiltonian can be mapped to a spinless fermion system by Jordan-Wignertransformation:
∏
∏
1
Si+ =
(1 − 2c†j cj )c†i , Si− =(2.3)(1 − 2c†j cj )ci , Siz = c†i ci − ,
2
0≤j<i
0≤j<i

Hfermion = −

N
−1 (
∑
i=0

(
)(
))
)
1( †
1
1
†
†
†
c ci+1 + ci+1 ci + ∆ ci ci −ci+1 ci+1 −
.
2 i
2
2

(2.4)

Boundary condition for spinless fermions is determined by that for the XXZ chain, although the relationship between them is slightly complicated [16]. Here we ignore it9

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 2. Definition of the problem: Flux quench in the S = 1/2 XXZ chain

Figure 2.1: Phase diagram of the ground state of the XXZ chain. For all anisotropies ∆,Hamiltonian is integrable and the exact ground state energy is known.

because in numerical calculations we treat the thermodynamic limit of the system, where
the boundary condition does not matter. In the spinless fermion formulation, ∆ is thestrength of nearest density-density interactions.
At ∆ = 0, Hfermion becomes tight-binding Hamiltonian and can be solved by Fouriertransformation:
N −1)∑
1 ∑( †
Hfermion,free = −cos(k)c̃†k c̃k ,(2.5)
ci ci+1 + c†i+1 ci = −2 i=0k
1 ∑ −ikrc̃k := √
cr e
.
N r=0
N −1

(2.6)

The ground state is thus the Fermi sea of spinless fermions, Πk≤kF c̃†k c̃k |vacuum⟩ (kF is the
Fermi momentum). Furthermore, it is known [11] that the XXZ Hamiltonian is exactly
solvable by Bethe ansatz for arbitrary anisotropies ∆. According to the exact solutions,
Hamiltonian is gapless for −1 ≤ ∆ ≤ 1, gapped for ∆ < −1 (anti-ferromagnetic phase)and 1 < ∆ (ferromagnetic phase). Phase diagramof the ground state is shown in Fig. 2.1.
∑
In gapless phase, there is no magnetization i ⟨Siz ⟩ and correlation functions ⟨S0x Srx ⟩ and⟨S0z Srz ⟩ show power-law decay. In addition, the model exhibits TLL universality [9]. The
renormalized Fermi velocity vF and the Luttinger parameter K are exactly known,√ππ 1 − ∆2
, K=
.
(2.7)
vF =
2 arccos(∆)
2 arccos(∆)∑
In gapped anti-ferromagnetic phase, non-zero staggered magnetization i (−1)i ⟨Siz ⟩ appears and the magnetization is zero. All correlation functions decay exponentially inspace. The ground state is similar to classical Néel state |↑↓↑↓ · · ·⟩. However, the true
ground state energy is lower than that of the Néel state and neighboring correlation−⟨Si+ Si+1⟩ is finite. Only in the limit of ∆ → −∞, the ground state becomes the Néelstate. In gapped ferromagnetic phase, the ground state is simply |↑↑ · · ·⟩ or |↓↓ · · ·⟩ forall ∆ > 1.
An important feature of the XXZ chain is integrability. The S = 1/2 XXZ chain is
integrable for any ∆ and its infinitely many conserved quantities {Qi }i can be calculatedby algebraic Bethe ansatz [11]. It is one of the non-trivial integrable models in the sensethat it is not equivalent to a model of free particles.
Previous studies on quantum quench in XXZ chain
The XXZ chain is one of the simplest models in one-dimensional quantum systems that
include interactions. For this reason, various types of quantum quench in the XXZ chain10

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 2. Definition of the problem: Flux quench in the S = 1/2 XXZ chain

Figure 2.2: Hamiltonian Hθ represents a ring of the XXZ chain pierced by magnetic flux
ϕ = θN .

have been studied. We review some of them here.
Barmettler et al. [17] studied quantum quench in the XXZ chain where initial stateis prepared as the Néel state |↑↓↑↓ · · ·⟩. They numerically calculated the dynamics ofstaggered magnetization by the iTEBD method. They found oscillatory dynamics in gapless phase (|∆| < 1) and nonoscillatory dynamics for large anti-ferromagnetic anisotropy(∆ ≪ −1). Furthermore, they showed that at the critical point (∆ = −1) the relaxation
time takes a minimum in contrary to the notion of critical slowing down, from which we
expect a maximum. Fagotti et al. [18] investigated quantum quench in the XXZ chain bytaking various initial states (Néel states, dimer states, etc.). They calculated the dynamicof some correlators such as ⟨S0x S1x ⟩ or ⟨S0z S2z ⟩ by time-dependent DMRG and iTEBD, andcompared the long-time limit of the correlators with analytical predictions from GGE.
They found the numerical data agree with the prediction from GGE well. The authors ofRef. [19] and Ref. [20] considered the so-called interaction quench in the XXZ chain. Theinteraction quench is a quantum quench of the anisotropy ∆ which is changed from zero
(free fermion point) to some finite value. They focused on the Z-factor of momentum
distribution of spinless fermions and kinetic energy of the system, and concluded thatthose quantities after the quench are well described by the eﬀective TLL. Their results
are in contrast with our present study, as we will see later. Finally, we comment thatanalytical methods based on the Bethe ansatz techniques to study quantum quench in
the XXZ chain are proposed by Ref. [21, 22].

2.1.1

XXZ chain with flux

The XXZ Hamiltonian with flux is defined as)∑ (11 −iθ − +iθ + −z z
e Si Si+1 + e Si Si+1 + ∆Si Si+1 .
Hθ = −
2
2i
In the spinless fermion representation,
(
)(
))
∑ (1
1 −iθ †
1
1
†
†
iθ †
Hθ = −
e ci ci+1 + e ci+1 ci + ∆ ci ci −ci+1 ci+1 −
.
2
2
2
2
i

(2.8)

(2.9)

θ can be viewed as Aharonov-Bohm phase on hopping of particles, or vector potential.
Hamiltonian Hθ corresponds to a situation where a ring of the XXZ chain is pierced bymagnetic flux ϕ = θN (Fig. 2.2). This Hamiltonian has been studied in the context of
the Drude weight (DC-conductivity) of the XXZ chain [23], spectral flow of the energy
levels [24] and the Berry’s phase induced by adiabatic insertion of the flux [25].
11《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 2. Definition of the problem: Flux quench in the S = 1/2 XXZ chain
For a finite size system, the presence of magnetic flux can be converted to twisted
boundary condition by Lieb-Schultz-Mattis (LSM) operator [16]
(
)
N
−1
∑
Uθ = exp −iθ
rSrz .
(2.10)
r=0

We can show easilyUθ† Hθ Uθ = H0,twisted = −

N
−1 (
∑
i=0

)
1 − +
1 + −
±
z z
Si Si+1 + Si Si+1 + ∆Si Si+1 , SN
= S0± e∓iθN . (2.11)
22This is zero-flux XXZ Hamiltonian with twisted boundary condition. An exact solution
is known [26] for the twisted boundary condition case and the ground state energy in thethermodynamic limit is the same as that of the XXZ chain without twist.
If θN ≡ 0 mod 2π, there is no eﬀects of twist even after the unitary transformationby LSM operator. In this case Uθ† Hθ Uθ = H0 holds, and therefore the ground state of Hθ
is given by that of H0 :|θ⟩GS = Uθ |θ = 0⟩GS .
(2.12)
Note that the ground state of Hθ and H0 is diﬀerent from each other although the energy
spectrum of them are the same.
In the spinless fermion language, LSM operator is written as
(
)
N
−1
∑
†
Uθ = exp −iθrcr cr .
(2.13)
r=0

Here we ignored an unimportant constant factor in Uθ . Uθ acts on fermion operators as
a boost in momentum space,Uθ† c̃k Uθ = c̃k+θ , Uθ† c̃†k Uθ = c̃†k+θ .

(2.14)

We will use this property later.
In this thesis, we study a spin current of the S = 1/2 XXZ chain. The definition ofspin current Jˆ is)
1 ∑1 ( + −+Jˆ :=ji , ji =Si Si+1 − Si− Si+1,
(2.15)
N i
2i
which is taken to satisfy1 the continuity equationd z!
Si = [iH0 , Siz ] = ji − ji−1 .
dt(2.16)We define Jˆ as an intensive quantity because we consider the thermodynamic limit N →∞ in numerical calculations. In spinless fermion notation, Jˆ is expressed as)
1 ∑( †1 ∑Jˆ =ci ci+1 − c†i+1 ci =sin(k)c̃†k c̃k .
(2.17)
2iN i
N k
1This definition of the spin current is for Hamiltonian H0 (zero-flux case). In the presence of flux,
the continuity equation changes and we have to slightly modify the definition of the spin current if we−+want to make it hold, jiθ = (2i)−1 (eiθ Si+ Si+1
− e−iθ Si− Si+1).

12

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 2. Definition of the problem: Flux quench in the S = 1/2 XXZ chain

Figure 2.3: Quench protocol of the flux quench. In our study, we take θ < 0 so that initial
spin currents become positive.

This is exactly the current of spinless particles. The spin current is not a constant of
motion in the presence of interactions (∆ ̸= 0):
ˆ =− ∆[H0 , J]
2iN

∑(

+
−
+ Si− Si+1
Si+ Si+1)( z)z
Si−1 − Si+2.

(2.18)

i

Non-conservationof the spin current is due to the Umklapp scattering by the term
∑ z z
∆ i Si Si+1 .

2.2

Flux quench in the S = 1/2 XXZ chain

Let us define flux quench in the S = 1/2 XXZ chain, a main problem of this study. The
flux quench is a quantum quench of phase (or flux) θ in the XXZ chain,{)∑ (1θ (t < 0)1 −iθ(t) − +z ziθ(t) + −H(t) = −.
e Si Si+1 + e
Si Si+1 + ∆Si Si+1 , θ(t) =2
20 (t ≥ 0)i(2.19)The protocol of quench is drawn in Fig. 2.3. At first (t < 0), the system is taken to be
the ground state. Suddenly the flux is turned oﬀ at t = 0, and then the system starts
evolving in time. If we regard2 −θ as a vector potential for spinless fermions, changing
the vector potential in time is equivalent to imposing electric field to fermions. In our
quench protocol, the electric field is delta function,
E(t) = −d(−θ(t))= −θδ(t).
dt(2.20)

Therefore, at t = 0 the (spin) current is induced by the flux quench.
Equivalently, we can regard initial spin current at t = 0 as a persistent current in
the ground state of the XXZ chain with flux. This is manifested in the spinless fermion
picture. As we discussed in the last section, when initial flux ϕ = θN is zero in modulo
2π the ground state of pre-quench Hamiltonian is written as Uθ |θ = 0⟩GS . Recalling thatUθ acts on fermion operators as a boost in momentum space, we know that Uθ |θ = 0⟩GS
In a standard notation,∑ a vector potential A is defined as HA ≈ HA=0 − Aj (j is current operator).
In our case Hθ = H0 + θ i ji + O(θ2 ), and therefore we define the vector potential as −θ.
213

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 2. Definition of the problem: Flux quench in the S = 1/2 XXZ chain
1
0.9
0.8

θ= 0
θ= -π/6
θ= -π/3

0.7

nq

0.6
0.5
0.4
0.3
0.2
0.1
0
-1 -0.8 -0.6 -0.4 -0.20

0.2 0.4 0.6 0.8

1

Momentum q /πFigure 2.4: Initial momentum distribution nq = ⟨c̃†q c̃q ⟩.

is the shifted Fermi sea by momentum −θ (initial momentum distribution nq = ⟨c̃†q c̃q ⟩ is
ˆ = N −1 ∑ sin(q)nq , the shiftedshown in Fig. 2.4). Since the (spin) current is given by ⟨J⟩qFermi sea has non-zero persistent current, which is the origin of initial current of the flux
quench.
In this study, we investigate the dynamics of spin current after the flux quench by
numerical calculations. We focus on two points:(i) What is the fate of spin current in the long-time limit? As we will see below, fromlinear response theory in θ the long-time limit of the current is given by the Drude
weight of the XXZ chain. How dose nonlinearity appear for large θ?
(ii) How is the dynamics after the quench? In gapless phase there is no a priori energy(or time) scales. What determines the time scale of the dynamics?
The main results of our study are presented in chapter 4. In the rest part of this section,
we discuss possibilities of experimental realizations of the flux quench and the analysis bylinear response theory on the flux quench.
Experimental realizations of the flux quench
Before going to analyze the flux quench by linear response theory, we comment on possible
experimental realizations of the flux quench. The flux quench can be viewed as a shiftof momentum at t = 0 in the spinless fermion picture. Therefore, it is equivalent to asituation where a moving lattice (its velocity gives an initial shift of momentum) stops
suddenly at t = 0. This quench process was experimentally realized [27] in bosonic
optical lattice. As for bosonic optical lattice, a quantum quench of artificial gauge field
is also proposed [28]. This is a more direct realization of the flux quench (although
it is a bosonic system). We expect that these kinds of quench will be performed forfermionic optical lattice in the near future. Furthermore, recently a optical lattice system
which is equivalent to the Heisenberg chain (∆ = 1) was realized and an experiment14

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 2. Definition of the problem: Flux quench in the S = 1/2 XXZ chain
on nonequilibrium transports in the chain was done [29]. Considering these substantial
experimental developments, we are optimistic about realizing our flux quench in opticallattice systems. We hope it provides us information that cannot be obtained by thenumerical approach employed in this thesis.

2.2.1

Analysis by linear response theory and the Drude weight

When θ is small, a term eiθ can be expanded as eiθ = 1 + iθ + O(θ2 ). Hamiltonian for theflux quench problem can be written as
H(t) ≈ Hθ − θJN Θ(t),(2.21)

where Θ(t) is step function. We can employ linear response (LR) theory in θ and obtain
ˆJ(t) := ⟨J(t)⟩,∫ ∞1J(t) =dωF (ω)σ(ω)e−iωt + O(θ2 ),(2.22)
2π −∞where F (ω) = −θ is Fourier transformation of the imposed electric field E(t) = −θδ(t)and σ(ω) is spin conductivity3
∫
N ∞
ˆ J]⟩
ˆ GS,0 .
σ(ω) =
dt eiωt ⟨[J(t),(2.23)
ω 0Here ⟨. . .⟩GS,0 means the expectation value for the ground state of H0 . The spin conductivity has anomalous zero-frequency component and it is called the Drude weight:
σ(ω) = 2πDδ(ω) + σreg (ω).

(2.24)

Thus the long-time limit of the current is given by D (we take θ < 0 in this study),J(t = ∞) = −Dθ = D|θ|.

(2.25)The Drude weight D is related to change of the ground state energy as inserting small
flux Φ into the system, and was calculated exactly by the Bethe ansatz techniques [23],D=

π sin µ, µ = arccos(−∆).
4 µ(π − µ)(2.26)

D is non-zero only for −1 ≤ ∆ < 1. In gapped phase (|∆| > 1) and ∆ = 1, D is zero.
Hence LR theory predicts that J(t = ∞) is non-zero only for −1 ≤ ∆ < 1.

2.3

Previous studies

In this section, we review previous studies related to our flux quench, although motivations
and setups of those studies are diﬀerent from ours. We introduce two kinds of studies
related to our problem; One is about the long-time limit of the current. The other isabout a bosonic version of flux quench in the context of dynamical stability of superfluid.
3

Non-perturbative Hamiltonian for LR theory in this case is Hθ , not H0 . Therefore expectation
values in LR formula should be taken for the ground state of Hθ . However, in O(θ) level we can obtainthe same results if expectation values are taken for the ground state of H0 .

15

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 2. Definition of the problem: Flux quench in the S = 1/2 XXZ chain

Figure 2.5: Left panel: The spin current after the flux quench for θ = −π/2, ∆ = −0.5
(red line). Note that in their notation the time scale and definition of the current is
twice as ours. They used time-dependent Lanczos method for L = 26 system [6]. Rightpanel: Estimation of the Drude weight by numerical data J(t) (t∗ /t0 is the time at which
estimation was done and ϕ0 is initial fiux). Blue line shows exact values of the Drude
weight. Note that their definition of ∆ is diﬀerent from ours by sign. ∆ > 0 in the figure
corresponds to ∆ < 0 in our notation. They used time-dependent DMRG for L = 100
system [28].

2.3.1Previous studies on the long-time limit of the current

Mierzejewski et al. [6] studied the flux quench in spinless fermion formulation, exactly thesame as our setup, to illustrate the breakdown of GGE. They showed that steady states
after the flux quench cannot be described by GGE. This result is the first example thatclearly denies GGE, and quite seminal for understanding the thermalization phenomena
in integrable systems. Their argument
was simple; conserved quantities Qi contained∑in GGE (ρGGE ∝ exp(−βH − i λi Qi )) for the XXZ chain have even parity under theparticle-hole transformation (ci → (−1)i c†i ). The spin current Jˆ has odd parity under
the same transformation. Therefore the expectation value of the spin current for any
ˆ = 0. Nevertheless, as LR theory predicts and numerically
GGE must be zero, Tr(ρGGE J)
verified by the authors, the long-time limit J(t = ∞) is obviously non-zero in some cases(left panel of Fig. 2.5). This is a clear illustration of the breakdown of GGE. The steadystates after the flux quench cannot be expressed by GGE as well as Gibbs (thermal)
states. From the viewpoint of our motivations for the flux quench, we comment that they
did not study J(t = ∞) systematically varying θ and ∆.
Moreover, a recent publication [28] mainly studying the flux quench in the BoseHubbard model mentioned the long-time limit of the current in the XXZ chain. They
dealt with the flux quench in the XXZ chain as a special limit of the Bose-Hubbard model.
They plotted (right panel of Fig. 2.5) estimation of the Drude weight from numerical dataJ(t) at late times (we do not show the details here). If LR theory is exact, all points inthe plot must be on the line. The deviation from the line (LR theory) is large for ∆ > 0(in our notation, ∆ < 0). They argued that this tendency is due to a finite-size eﬀect of
16

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 2. Definition of the problem: Flux quench in the S = 1/2 XXZ chain

Figure 2.6: The spin current after the flux quench in one-dimensional Bose-Hubbard
model [30]. They used iTEBD.
simulations and no further discussion was made. As we will see in section 4.2, we confirm
this tendency but claim that the origin is not a finite-size eﬀect; it is an intrinsic behaviorof the system relating to superfluidity of the system.

2.3.2

Bosonic flux quench

Dynamics of the Bose-Hubbard model has been studied extensively because of its direct
experimental realizations by optical lattices. In this context, the stability of currents insuperfluid regime of the Bose-Hubbard model has been theoretically and experimentallyinvestigated. The authors of Ref. [30] considered a problem of abrupt momentum boost(equivalent to quench of flux) for the ground state of one-dimensional Bose-Hubbard
model,) U∑∑( ††bi bi+1 + bi+1 bi +HBH = −Jni (ni − 1).
(2.27)
2 ii
bi is an annihilation operator of boson on site i, ni = b†i bi is the number of boson on site
i. Definition of the current operator is the same as in our spinless fermion case,
)
J( †
†
ji =
b bi − bi bi+1 .
i i+1(2.28)

They took initial state as the boosted ground state of HBH by momentum ka. The
momentum ka corresponds to θ in our flux quench problem. The dynamics they obtained
is shown in Fig. 2.6. For small momentum ka and small interaction u = U/J, the current
decays slowly and sometimes remains non-zero in the long-time limit. The behaviors of
the current are similar to our results, although the model under investigation has diﬀerent
particle statistics (boson and spinless fermion). Besides, there are no oscillations in the
dynamics of bosonic current, in contrast to our results. We will see other results of bosonicflux quench in section 4.4 to discuss the relaxation time of the current.

17

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3
Numerical Methods
We employ infinite Time-Evolving Block Decimation (iTEBD) method for numerical calculations in this study. In this chapter, we review numerical methods. The iTEBD isbased on matrix product states (MPS), a kind of tensor networks. An overview of tensor
networks is given in the first section. In the second section, background of tensor networks
is reviewed. For a comprehensive and pedagogical review of tensor networks, see Ref. [31].
In the third section, the basics of matrix product states with some detailed proofs arepresented. In the last two sections, Time-Evolving Block Decimation (TEBD) method
and iTEBD method are explained. This chapter is almost independent of other contents
in this thesis, so readers who are not interested in details of numerical methods can skip
here.

3.1

Introduction

Quantum many-body systems are at the heart of studies in condensed matter physics.
Interactions and correlations between a single degree of freedom (particles, spins, etc.)
evoke interesting phenomena which cannot be attributed to the property of such a single
constituent: for example, superconductivity, fractional quantum Hall eﬀect, non-Fermi
liquids and quantum spin liquids. In order to study these phenomena, a simple model
is often proposed, such as Hubbard model for superconductivity or Heisenberg modelfor quantum magnets and spin liquids. However, it is far diﬃcult to solve these models
exactly, except for very limited examples. Therefore one often resorts to take numericalapproaches to see the nature of these models. Nevertheless, numerical calculations have
limitations, too. Exponential growth in a dimension of Hilbert space with a system size
restricts available numerical calculations for large systems. For example, the dimensionof Hilbert space of N site lattice model containing M spinless fermions is
( )
N!
N
,
(3.1)
dim H ==M
M !(N − M )!
which is exponentially large in N for a fixed particle density M/N . There are many analytical and numerical methods to overcome this diﬃculty: mean-field theory, perturbative
expansions and diagrammatic techniques, renormalization group, dynamical mean-field
theory, quantum Monte Carlo, and so on. All of these methods have advantages and18

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
disadvantages. For example, quantum Monte Carlo can deal with larger system size thanexact diagonalization can, although the so-called negative sign problem [32] prevents us
from applying the method to frustrated quantum magnets and fermionic systems.
Tensor network methods are one kind of approaches to tackle the problem of large
dimension of Hilbert space. They do not have the negative sign problem and can beapplied to bosonic, fermionic and spin systems in an almost same way1 . One importantconcept in tensor networks is entanglement. If given quantum state has a small amount ofentanglement, tensor networks can express it by a smaller amount of information comparedwith a usual expression using wave function. Area law of entanglement entropy assuresthat the ground states of Hamiltonian which we often deal with have smaller amounts ofentanglement than the vast majority of states in Hilbert space. Tensor network methods
have obvious criteria for their applicability, and the error in computation can be estimatedin a clear way (especially for one-dimensional systems).
Famous examples of tensor network methods are density matrix renormalization group(DMRG) [14], time-evolving block decimation (TEBD) [15], projected entangled pairstates (PEPS) [34] and multi-scale entanglement renormalization ansatz [35]. DMRG,
proposed by S. White in 1992 [14], is one of the most successful numerical methods for
one-dimensional quantum systems. Although DMRG was proposed independently fromthe notion of tensor networks, the close relationship between them was realized later (fora review, see Ref. [36]).
Finally, it should be stressed that tensor networks are not only numerical approaches
but also conceptually new tools to understand quantum phases. Recently, quantum
phases which cannot be classified in terms of Landau’s symmetry breaking were noticed and named topological ordered phases [37, 38]. Symmetry protected topologicalphases, one variant of topological ordered phases, can be classified by MPS and PEPS
(one-dimensional and two-dimensional version of tensor networks) [39–41].

3.2

Tensor networks

In tensor network methods, a wave function of quantum state is regarded as a large-rank(many legs) tensor and expressed as contracted sum of many small-rank tensors. One of
the advantages of the tensor network methods is eﬃciency for expressing quantum states.
In this section, we first introduce the notation and diagrammatic representation of tensor
networks by taking a toy example of tensor network. Next, the area law of entanglement
entropy, which is the most important fact behind tensor networks, is explained. Finally,
the relationship between tensor network states and area-law states is discussed.

3.2.1

Basics of tensor networks

Consider rank-4 tensor Aαβγδ as a simple toy example of tensor networks. Indices α . . . δ
run from 1 to χ. χ is called the bond dimension2 and can be diﬀerent for each index of the1

Tensor networks for fermionic systems have a problem of sign originating from order of fermion
creation operators. However it is not so serious as the negative-sign problem in quantum Monte Carlo [33].
2
In DMRG literature, the bond dimension is written as m or D. In TEBD literature, χ is usually
used.

19

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

Figure 3.1: Diagrammatic representation of Eqn. (3.2).

tensor (for simplicity here we make χ the same among indices). Assume that the tensor
A is decomposed into two rank-3 tensors B and C,Aαβγδ =

χ̃∑

Bαβϵ Cϵγδ .

(3.2)

ϵ=1

We introduce a diagrammatic representation for calculations in tensor networks, because
it is helpful to understand complicated equations with many indices. Equation (3.2) isdiagrammatically written as Fig. 3.1. Tensors are expressed in diagrams such as circles
and the number of legs shows the rank of the tensor.
Let us consider this decomposition from the viewpoint of an amount of information.
Aαβγδ has χ4 components, while Bαβϵ and Cϵγδ have χ̃χ2 components. If new bonddimension χ̃ is equal3 to χ, the number of components to describe the right hand sideof Eqn. (3.2) is χ3 + χ3 = 2χ3 . Hence when χ is large, the expression in the right handside of Eqn. (3.2) is more eﬃcient than the left hand side (raw tensor) in terms of theamount of information required to store the tensors. This is an illustration of one of theadvantages of tensor networks.
In tensor network methods, we decompose a wave function of a quantum state (it isa rank-L tensor when system size is L) into many small-rank tensors as in the previous
toy example. That decomposition can greatly reduce the required memory storage andcomputational cost in numerical calculations. When attempting to use decompositions
like the previous example for general wave functions in Hilbert space, one may pose thesequestions:(a) how to calculate decompositions for general wave functions and obtain tensor network descriptions of them?
(b) in which conditions the bond dimension of small-rank tensors decomposing a wave
function is small so as to make the tensor network description eﬃcient?
For one-dimensional systems, we will answer both questions in section 3.3. For generaldimensions we discuss only the latter question in the rest part of this section. The answerof the latter question is entanglement of a quantum state to be represented in tensor
networks.
3

In general, χ̃ becomes larger than χ. As we will see later, the condition in which we can obtain smallbond dimension χ̃ is related to entanglement of a tensor to be decomposed.

20

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

3.2.2

Entanglement entropy and area law

As we will see in section 3.3, the bond dimension χ of tensor network methods can be
estimated by entanglement entropy of a state to be expressed. Entanglement entropy is
one of the most famous measures of entanglement in quantum systems. In this subsection,
we show the definition of entanglement entropy and introduce the area law of entanglement
entropy.
Consider a system whose Hilbert space H can be written as a tensor product of Hilbert
spaces of two subsystems A and Ā, H = HA ⊗HĀ . For lattice systems, it is always possible
to do this decomposition. Entanglement entropy SA of a pure state |ψ⟩ ∈ H is defined as
SA := −TrA (ρA ln ρA ) ,(3.3)

where ρA = TrĀ (|ψ⟩ ⟨ψ|) is the reduced density matrix of subsystem A. Entanglement
entropy is von-Neumann entropy of the reduced density matrix ρA ; this means that theentanglement entropy SA quantifies information contained in ρA . If two subsystems A
and Ā are entangled, tracing out the subsystem Ā reduces the information in the subsystem A and then SA becomes large. This is why entanglement entropy is a measure ofentanglement.
Next we introduce the area law of entanglement entropy. This law assures that tensor
networks can be useful to represent the ground states of Hamiltonians that we often
encounter.
Area law of entanglement entropy (review: [42]) Consider a lattice system. Groundstates of gapped, local Hamiltonians satisfy
SA ∼ const. |∂A| + o(|∂A|),(3.4)

where SA is entanglement entropy for subsystem A and |∂A| is the area of boundary ofA, ∂A. By local we mean that supports of each term in the Hamiltonian are finite, notscaling as system size L. Furthermore, in one-dimensional systems, the ground states ofgapless and local Hamiltonians also exhibit area law with logarithmic corrections,SA ∼ const. |∂A| ln(|A|) + O(|∂A|).

(3.5)

A physical interpretation of the area law is following. The correlation length ξ of theground state of gapped and local Hamiltonian is finite [43]. Then each site of the systemcannot be entangled over the distance ξ, which means that entanglement between two
subsystems A and Ā is mostly carried by the interface of them (Fig. 3.2). Therefore, SA
is roughly estimated to be proportional to |∂A|. There are rigorous proofs for the area law
in one-dimensional quantum systems under some conditions [44]. For higher-dimensional
systems, rigorous results are limited.
Since we mainly consider one-dimensional systems in this thesis, let us explain thearea law in one dimension more concretely. Consider a L site lattice system (chain) in
one dimension and a gapped and local Hamiltonian on it. The area law implies thatentanglement entropy of the ground state for a half chain (one half of the system) does
not scale as system size L,
Shalf chain ∼ O(1),
(3.6)
21

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

Figure 3.2: Division of a lattice system into two subsystems A and B = Ā. If thecorrelation length ξ is finite, only the sites close to the boundary of A within the distanceξ (shaded area) can contribute to the entanglement between A and Ā. Therefore theentanglement entropy SA is roughly proportional to the size of boundary |∂A|. Figure istaken from Ref. [45].

even though the size of the half chain scales as L/2. As we will see in section 3.3, this result
greatly aﬀects applicability of tensor networks in one dimension. For gapless Hamiltonian,
the area law indicates Shalf chain ∼ ln(L). Moreover, we can exploit conformal field theory
to calculate entanglement entropy for this case; the result is [46]Shalf chain ∼

c
ln L,
6

(3.7)

where c is the central charge of conformal field theory describing the critical (gapless)ground state.

3.2.3

Ground states are well expressed by tensor networks

In addition to the area law, it is known [47] that if we randomly choose a state in Hilbert
space the state has entanglement entropy which obeys volume law, SA ∼ |A|, with probability one. In other words, almost all of the states in Hilbert space obey the volume lawof entanglement entropy. Therefore, the ground state of gapped and local Hamiltonian,which obeys the area law, sits on the very tiny corner of the Hilbert space (Fig. 3.3).
Furthermore, it is proved [48] that the amount of quantum many-body states that canbe generated by arbitrary time-dependent local Hamiltonians in a time that scales polynomially in the system size is exponentially small in Hilbert space. The vast majority ofstates in Hilbert space is not physical or reachable in this sense. Considering these twofacts, we can say that the whole Hilbert space is too large to describe the physics aroundthe ground state.
Tensor network states with constant (not scaling as system size) bond dimensions
exhibit the area law (although we do not show a explanation [31] here). Hence they aregood candidates to describe states around the ground states. Actually, it is proved [45,49]
that tensor network states with constant bond dimensions can express states around theground states eﬃciently. This is an important conceptual point of tensor network methods.
22《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

Figure 3.3: Almost all of the states in Hilbert space obey the volume law. States which
exhibit area law (including the ground state) occupy a tiny amount of Hilbert space.
Tensor network states can describe such area-law states eﬃciently.

3.3

Matrix product states (MPS)MPS is a kind of tensor networks for one-dimensional quantum systems. MPS is a basisof two successful numerical methods for one-dimensional quantum systems, DMRG and(i)TEBD. In this section, we review MPS. First we define MPS and introduce Schmidt
value decomposition, which is important to understand the structure of MPS. Next weshow that any pure state in a finite lattice system can be expressed in a MPS form. Then
canonical conditions of MPS are stated and calculations of physical quantities in MPS
language are described. We show examples of MPS in the last part of this section.
Let us consider a lattice system of L site whose Hilbert space consists of L times( )⊗L
tensor product of local Hilbert space, H = Cd ⊗ Cd ⊗ · · · ⊗ Cd = Cd
. The basis ofdlocal Hilbert space is written as {|i⟩}i=1 and d is called the local dimension. For spin-S
systems, {|i⟩} = {|sz = −S⟩ , . . . , |sz = S⟩} and d = 2S + 1. For spinless fermion systems,{|i⟩} = {|0⟩ , |1⟩} and d = 2. For bosonic systems, {|i⟩} = {|0⟩ , |1⟩ , · · · } and d = ∞ (weset the cut-oﬀ dmax on d in actual numerical calculations). Given pure state |ψ⟩ in thissystem can be expressed by a wave function ci1 i2 ···iL ,
d
∑

|ψ⟩ =

ci1 i2 ···iL |i1 i2 · · · iL ⟩ .

(3.8)

i1 i2 ···iL =1

Here, |i1 i2 · · · iL ⟩ is a tensor product of the basis of local Hilbert space and constitutes
the complete orthonormal basis of H. Matrix product state {Γ, λ} for |ψ⟩ is defined4 as
|ψ⟩ =

d
∑

χi
∑

[2]
[L−1]iL−1 [L−1] [L]iL
1 [1] [2]i2Γ[1]iα1 λα1 Γα1 α2 λα2 · · · ΓαL−2 αL−1 λαL−1 ΓαL−1 |i1 i2 · · · iL ⟩ .

(3.9)

i1 i2 ···iL α1 ···αL−1 =1
4

There are two major notations of MPS. We used Γλ-notation here. The other notation is A-notation
|ψ⟩ =

d
∑

Di
∑

[2]i2
[L−1]iL−1 [L]iL1
A[1]iα1 Aα1 α2 · · · AαL−2 αL−1 AαL−1 |i1 i2 · · · iL ⟩ ,i1 i2 ···iL α1 ···αL−1

which is often used in the DMRG literature.

23

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

Figure 3.4: A diagrammatic representation of MPS.

This seems complicated at first sight, but a diagrammatic representation in Fig. 3.4 ishelpful to understand the structure of tensors. Note that the bond dimension χi depends
[l]il
[l]
on site i. If we regard tensors Γαl−1
αl as χl−1 × χl matrices and λαl as χl × χl diagonal
[1]i
[L]iL
matrices for each il (Γα1 1 and ΓαL−1
as 1 × χ1 (L−1) vectors), |ψ⟩ can be written as thesum of products of matrices,
|ψ⟩ =

)T
∑ (
⃗Γ[1]i1 λ̂[1] Γ̂[2]i2 λ̂[2] · · · Γ̂[L−1]iL−1 λ̂[L−1]⃗Γ[L]iL |i1 i2 · · · iL ⟩ .

(3.10)

i1 i2 ···iL

This is the origin of the name Matrix Product States.
One of the advantages of using MPS representations for quantum states is that MPS
requires the smaller number of parameters than the usual wave function descriptions
ci1 i2 ···iL do. In order to retain the information of wave function ci1 i2 ···iL , dL ∼ exp(L)parameters are needed. On the other hand, at most dχ2 L + Lχ parameters are needed
to store the tensors Γ and λ of MPS (χ := maxi (χi )). The value χ determines eﬃciencyof MPS. As we will discuss later, χ is roughly estimated by the entanglement entropyof |ψ⟩. For general states in Hilbert space, or volume-law states, χ scales as exp(L) andtwo representations (wave functions and MPS) are equivalent in terms of the amountof information required to express the states. However, for the ground state of local
Hamiltonian in one dimension, or area-law states, χ does not scale as L (in gappedsystems) or scales as polynomials of L (in gapless systems). In both cases, we can bounddχ2 L+Lχ ≲ poly(L). This is far smaller than exp(L), and therefore MPS is more eﬃcientthan the wave functions to express the ground state of local Hamiltonian.

3.3.1

Schmidt decomposition

Schmidt decomposition of a state in Hilbert space is essential to construct (canonical) MPS
representations for general quantum states. In order to define Schmidt decomposition, we
first introduce singular value decomposition for general matrices [50, §2.5].
Singular Value Decomposition (SVD). Let M be a m1 × m2 matrix. Singular value
decomposition of M is defined as
M = U SV,
(3.11)
where U and V are m1 × min(m1 , m2 ) and min(m1 , m2 ) × m2 matrices which satisfy
U † U = 1 and V V † = 1, S is min(m1 , m2 )×min(m1 , m2 ) diagonal matrix with non-negative
diagonal elements Sα ≥ 0. S is unique up to permutation of the diagonal elements. When24

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
m1 = m2 , U U † = 1 and V † V = 1 also hold and U and V are unitary matrices.
Diagonal elements of S are called singular values of matrix M . We can define {Sα } indescending order and omit zero(s)5 from them, S1 ≥ S2 ≥ · · · ≥ Sχ > 0. χ depends onmatrix M and satisfies 1 ≤ χ ≤ min(m1 , m2 ). Then we can write down SVD as
Mij =

χ
∑

Uiα Sα Vαj , Sα > 0.

(3.12)

α=1

We will use this notation of SVD in the following.
A
By using SVD, we define Schmidt decomposition of a state |ψ⟩. Let {|i⟩A }di=1
anddB{|j⟩B }j=1 be complete orthonormal bases of two subsystems A and B, where A ∪ B isequal to a total system. |ψ⟩ can be expanded as∑

dA ,dB|ψ⟩ =cij |i⟩A |j⟩B .

(3.13)

i,j=1

Coeﬃcient matrix cij is decomposed by SVD as
cij =

χ
∑

Uiα λα Vαj , λα > 0.

(3.14)

α=1

Here we write the singular values as λ. Putting this decomposition into Eqn. (3.13), one
obtainsχ∑ ∑dA ,dB|ψ⟩ =Uiα λα Vαj |i⟩A |j⟩B =i,j=1 α=1

where
|α⟩A :=

λα |α⟩A |α⟩B ,(3.15)

α=1dA∑

Uiα |i⟩A , |α⟩B :=

i=1

∑χ

χ
∑

dB
∑

Vαj |j⟩B .

(3.16)j=1|ψ⟩ = α=1 λα |α⟩A |α⟩B , {|α⟩A,B } and {λα } are called Schmidt decomposition, Schmidt†
vectors and Schmidt values, respectively. Because of the conditions U∑
U = 1 and V V † =χ
1, {|α⟩A,B }α=1 are orthonormal (⟨α|β⟩A,B = δαβ ). From this property, α (λα )2 = ⟨ψ|ψ⟩ =1 follows if the state |ψ⟩ is normalized. Besides, the reduced density matrix ρA is written
down in a simple form,()∑∑
ρA = TrB (|ψ⟩ ⟨ψ|) = TrBλα λ∗β |α⟩A |α⟩B ⟨β|A ⟨β|B =(λα )2 |α⟩A ⟨α|A . (3.17)
αα,βEntanglement entropy SA is
SA = −TrA (ρA ln ρA ) = −∑(λα )2 ln(λα )2 .
α5Here we assume M ̸= 0.

25(3.18)

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
Note that Schmidt value decomposition depends on the way of dividing total Hilbert
space.
The number of Schmidt values, χ, reflects the quantum nature of |ψ⟩. When χ = 1,
|ψ⟩ is a direct product of two states in each subsystem and can be considered as classical.
χ > 1 means that two subsystems A and B are entangled and that |ψ⟩ cannot be written
as a classical direct product of two states. If two subsystems are maximally entangled, χ =√min(dA , dB ), λα ≡ 1/ χ and SA = ln(χ). Conversely, when we fix χ, the states written
in the form of Schmidt decomposition (3.15) can represent states that have entanglement
entropy SA = ln χ at most. This argument implies a rough relation6 between χ and SA ,
χ ∼ eSA .

3.3.2(3.19)

Construction of MPS for arbitrary states in a system

To prove that any quantum state can be expressed in a MPS form, we repeatedly apply
SVD to a wave function of a state. Given quantum state |ψ⟩ can be expanded as∑|ψ⟩ =ci1 i2 ···iL |i1 i2 · · · iL ⟩ .
(3.20)
i1 i2 ···iL

As a first step, we regard i2 · · · iL as one index and perform SVD for matrix ci1 ,(i2 ···iL ) :
χ1
∑

SVD

ci1 ,(i2 ···iL ) =[1][1]
Ui1 α1 λ[1]α1 Vα1 (i2 ···iL ) =:

α1 =1

χ1
∑

[1]

1 [1]
Γ[1]i
α1 λα1 Vα1 (i2 ···iL ) .

(3.21)

α1 =1This is actually Schmidt decomposition of |ψ⟩ for the subsystem of site 1 and the subsystem of site [2 · · · L],
∑
∑
[1]
1 [1]
ci1 i2 ···iL |i1 i2 · · · iL ⟩ =
Γ[1]i
|ψ⟩ =
α1 λα1 Vα1 (i2 ···iL ) |i1 ⟩ |i2 · · · iL ⟩i1 i2 ···iL=
where|τα[1]1 ⟩ :=∑i1 ···iL α1χ1∑[1][2···L]λ[1]⟩,α1 |τα1 ⟩ |τα1
α1 =1[2···L]
1
Γ[1]i
⟩ :=
α1 |i1 ⟩ , |τα1i2 ···iL

i1
[1]

∑

[2···L]

Schmidt vectors {|τα1 ⟩}α1 and {|τα1and H[2···L] . Next we apply SVD to[1]

Vα1 (i2 ···iL ) |i2 · · · iL ⟩ .

(3.22)(3.23)⟩}α1 are orthonormal vectors in Hilbert spaces H[1]

[1]

[1]

V α1 i2 ···iL := λ[1]α1 Vα1 (i2 ···iL )(3.24)[1]by regarding α1 i2 as one index, V (α1 i2 )(i3 ···iL ) . We obtain Γ[2] and λ[2] by
χ2
χ2
∑
∑
[1]
SVD
[2]
[2] [2][2]i2 [2] [2]V (α1 i2 )(i3 ···iL ) =λ[1]U(α1 i2 )α2 λα2 Vα2 (i3 ···iL ) =:α1 Γα1 α2 λα2 Vα2 (i3 ···iL ) .
α2 =1
α2 =1
6

(3.25)

This is a physical, not mathematically rigorous relation. What we can say rigorously is that thevalue of entanglement entropy sets restrictions on the distribution of the singular values {λi }χi=1 , not on
the value χ itself. The detailed discussion is found in Ref. [44].

26《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
[2]i

[2]

[1]

Note that we define Γα1 α22 as U(α1 i2 )α2 /λα1 , which is always possible because of the property[1]of SVD, λα1 > 0. Using Γ[2] , we can define orthonormal vectors in Hilbert space H[12] andH[3···L] ,∑|τα[12]
⟩ :=
2

[2]i2
[1]
λ[1]
α1 Γα1 α2 |τα1 ⟩ |i2 ⟩ =i 2 α1|τα[3···L]
⟩
2

∑

:=

i3 ···iL∑

1 [1] [2]i2Γ[1]iα1 λα1 Γα1 α2 |i1 i2 ⟩ ,(3.26)

i1 i2 α1[2]Vα2 (i3 ···iL ) |i3 · · · iL ⟩ .

[12](3.27)[1]

Orthonormality of {|τα2 ⟩}α2 follows from that of {|τα1 ⟩}α1 and the property of matrix(U [2] )† U [2] = 1. Again, this is Schmidt decomposition of |ψ⟩ for the subsystems [12] and[3 · · · L],
χ2
∑
[12]
[3···L]
|ψ⟩ =
λ[2]
⟩.
(3.28)
α2 |τα2 ⟩ |τα2
α2 =1

One can prove this equation from the definitions of Γ[1] , Γ[2] , λ[1] and λ[2] .
We repeat this procedure of SVD until the end of system is reached (k = 2, · · · , L−2).
[k]

V (αk ik+1 )(ik+2 ···iL )

:=

[k]

λ[k]
αk V(αk ik+1 )(ik+2 ···iL )
∑

χk+1
SVD

=

[k+1]

[k+1]

U(αk ik+1 )αk+1 λ[k+1]αk+1 Vαk+1 (ik+2 ···iL )αk+1 =1

=:

∑

[k+1]

[k+1]ik+1 [k+1]λ[k]αk Γαk αk+1 λαk+1 Vαk+1 (ik+2 ···iL ) ,(3.29)αk+1⟩|τα[1···k+1]
k+1

:=

∑

[k+1]ik+1[1···k]⟩ |ik+1 ⟩λ[k]αk Γαk αk+1 |ταk

ik+1 αk∑=[k+1]ik+1
1 [1] [2]i2Γ[1]iα1 λα1 Γα1 α2 · · · Γαk αk+1 |i1 · · · ik+1 ⟩ ,(3.30)

i1 ···ik+1
α1 ···αk|τα[k+2···L]
⟩
k+1

:=

∑

ik+2 ···iL[k+1]

Vαk+1 (ik+2 ···iL ) |ik+2 · · · iL ⟩ .

(3.31)

Finally, we obtain a MPS representation of |ψ⟩
|ψ⟩ =

∑

ci1 ···iL |i1 · · · iL ⟩

i1 ···iL

=

∑

[1]

1 [1]
Γ[1]i
α1 λα1 Vα1 i2 ···iL |i1 · · · iL ⟩

i1 ···iL
α1

=

∑

[2]

[2]
1 [1] [2]i2Γ[1]iα1 λα1 Γα1 α2 λα2 Vα2 i3 ···iL |i1 · · · iL ⟩

i1 ···iL
α1 α2

=

∑

[L−1] [L]iL
1 [1] [2]i2Γ[1]iα1 λα1 Γα1 α2 · · · λαL−1 ΓαL−1 |i1 · · · iL ⟩

i1 ···iL
α1 ···αL−1

27(3.32)

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

Figure∑3.5: Diagrammatic representation of the construction of MPS for arbitary states|ψ⟩ = i1 ···iL ci1 ···iL |i1 · · · iL ⟩.
with Schmidt vectors for decomposing the system into H[1···k] ⊗ H[k+1···L]
|ψ⟩ =

χk
∑

[1···k]λ[k]⟩ |τα[k+1···L]⟩
αk |ταkk(3.33)αk =1|τα[1···k]
⟩
k

∑

=

[k]ik
1 [1]
Γ[1]i
α1 λα1 · · · Γαk−1 αk |i1 · · · ik ⟩ ,(3.34)

i1 ···ik+1
α1 ···αk−1∑⟩ =|τα[k+1···L]k[L]iL
k+1 [k+1]
Γ[k+1]iαk αk+1 λαk+1 · · · ΓαL−1 |ik+1 · · · iL ⟩ .

(3.35)

ik+1 ···iL
αk+1 ···αL−1

A diagrammatic representation of this construction of MPS is shown in Fig. 3.5.
Let us discuss what determines the bond dimension of MPS, χ := maxk (χk ). χkis the number of Schmidt values for dividing the system (into sites [1 · · · k] and sites)[k + 1 · · · L]. Therefore, from the property of SVD, χk ≤ min dim H[1···k] , dim H[k+1···L] =min(dk , dL−k ). This yields a general bound()χ = max χk = max min(dk , dL−k ) ≤ dL/2 .
(3.36)kk

Furthermore, the decompositions that we used to construct MPS are Schmidt decompositions indeed, so χ is related to the entanglement entropy of |ψ⟩. Recalling the relation (3.19), we know
S[1···k] ∼ ln χk ⇔ χk ∼ eS[1···k] .
(3.37)
Hence a physical estimation of χ is obtained:χ ∼ eSmax .

(3.38)

Here Smax is the maximum value of entanglement entropy of |ψ⟩ for all possible divisionsof the system.
28

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

Figure 3.6: Canonical conditions of MPS for 1 < k < L.

For a vast majority of states in Hilbert space, the volume law of entanglement entropy
Smax ∼ L holds and therefore χ scales as exp(L). As we already discussed, in that case
MPS representations require the same number of parameters to describe states as theusual wave functions do, and MPS is not eﬃcient to express states. However, for the
ground state of local Hamiltonian, the area law ensures Smax ∼ O(1) in gapped case andSmax ∼ (c/6) ln L in gapless case. This leads χ ∼ O(1) in gapped case and χ ∼ Lc/6 ingapless case. In both cases, the descriptions by MPS are far eﬃcient than that by wave
functions. This is why MPS is suitable and eﬃcient to describe the ground state of local
Hamiltonian.

3.3.3

Canonical condition for MPS

We have constructed MPS for any state in the previous subsection. However, there are
diﬀerent sets of matrices that express the same state. In other words, there is a choice of
freedom in MPS representations. Specifically, if MPS {Γ[k]ik , λ[k] }Lk=1 represents |ψ⟩, then()
e[k]ik = X [k−1] −1 Γ[k]ik X [k]Γ()
e[k] = X [k] −1 λ[k] X [k+1]λ(3.39)(3.40)is also a MPS representation7 of |ψ⟩ for arbitrary invertible χk × χk matrices X [k] .
Among many MPS representations for a given quantum state, we introduce canonicalconditions for MPS, which greatly reduce computational cost of physical quantities. The
canonical conditions of MPS are defined as8)∗(∑[k]ik [k]k [k]= δαβ ,
(3.41)λΓλ
Γ[k]iααk αkβαk αk∑ik ,αk)∗([k−1] [k]ik[k]ik= δαβ .
Γλ
Γ
λ[k−1]
αk−1 αk−1 βαk−1 αk−1 α(3.42)

ik ,αk−1

These two equations are called as right canonical condition and left canonical condition(diagrammatic representations are shown in Fig. 3.6). As we will see later, retaining
canonical conditions of MPS is crucial for eﬃcient computation of physical quantities.
()
e [1]i1 = Γ[1]i1 X [1] , Γe [L]iL = X [L−1] −1 Γ[L]iL .
For k = 1 and L, Γ
(
)
(
)∗
∗
∑
∑
[1]i
[1]i
[L]iL
[L]iL
8
For k = 1 and L, i1 Γα 1 Γβ 1 = δαβ ,ΓΓ= δαβ .
α
iL
β
7

29

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
It is proved that [51] MPS representations for given state that satisfy canonical conditions are unique up to permutation of the degenerate Schmidt values. The Γ and λobtained in the last subsection satisfy the canonical conditions (this is shown by straightforward calculations). Therefore, it is always possible to represent any quantum state ina canonical MPS.

30

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

3.3.4

Calculation of physical quantity using MPS

We will show MPS expressions of expectation values for one-site and two-site operators.
The canonical conditions of MPS are important to obtain simple formulae for the expectation values.

Expectation value for one-site operator[l]

Let Ô1 be a one-site operator acting nontrivially on site l and trivially on the other sites,[l][l]Ô1 = 1̂ ⊗ · · · ⊗ Ô1 ⊗ · · · ⊗ 1̂. The expectation value ⟨ψ|Ô1 |ψ⟩ can be written as[l]⟨ψ|Ô1 |ψ⟩)∗ ()∑ ( [1]j [1]
[L]j
[l]
[L]iL
1 [1]
=
Γβ1 1 λβ1 · · · ΓβL−1L
λ
·
·
·
Γ
Γ[1]i
α1
α1
αL−1 ⟨j1 · · · jL |Ô1 |i1 · · · iL ⟩
i1 ···iL ,j1 ···jL{α},{β}=∑ (()∗ ) (()∗ ) (
)∗ )
(
[l]
[1]i1
[1] [2]i2[l−1] [l]jl[1] [2]i2[l−1] [l]il
[l]
1
λ
Γ[1]i
Γ
λ
Γ
λ
Γ
·
·
·
λ
Γ
Γ
λ
⟨j
|
Ô
|i
⟩
λ
l
1 l
α1
α1 α1 α2αl−1 αl−1 αl αlβ1
β1 β1 β2βl−1 βl−1 βl βli1 ···iL ,jl{α},{β}=∑i2 ···iL−1 , jl{α},{β}

=

∑

)∗ )
(
)∗ ) ((([L−1]iL−1 [L−1]
[1]iL
[L]iL
[L−1]
L−1
Γ
Γ
λ
Γ
λ
× · · · Γ[L−1]i
αL−1
αL−2 αL−1 αL−1
βL−1
βL−2 βL−1 βL−1
(
)∗ )
(
(
)∗ ) ([l−1] [l]jl[l][1] [2]i2[1] [2]i2[l][l−1] [l]ilδα1 β1 λα1 Γα1 α2 λβ1 Γβ1 β2· · · λαl−1 Γαl−1 αl λαl ⟨jl |Ô1 |il ⟩ λβl−1 Γβl−1 βl λβl)∗ )(([L−1]iL−1 [L−1]
[L−1]
L−1
δαL−1 βL−1λΓλ× · · · Γ[L−1]i
αL−2 αL−1 αL−1
βL−2 βL−1 βL−1(()∗ )[l−1] [l]jl
[l]
[l]il
[l]
δα2 β2 · · · λ[l−1]
Γ
λ
⟨j
|
Ô
|i
⟩
λ
Γ
λ
· · · δαL−2 βL−2l1 l
αl−1 αl−1 αl αlβl−1 βl−1 βl βli3 ···iL−2 , jl
α2 ···αL−2
β2 ···βL−2

=···
=

∑

[l]il
[l]
λ[l−1]
αl−1 Γαl−1 αl λαl ⟨jl |Ô1 |il ⟩

(

[l]jl
[l]
λ[l−1]
αl−1 Γαl−1 αl λαl)∗.

(3.43)

il ,jl ,αl−1 ,αlA diagrammatic representation of this computation is shown in Fig. 3.7.
Therefore, we can calculate expectation values of one-site operators by using only threetensors of MPS, λ[l−1] , λ[l] and Γ[l] . The canonical conditions of MPS play an importantrole in the derivation of this formula.

Expectation value for two-site operator[l,l+1]Let Ô2
be a two-site operator acting nontrivially on site l and l + 1 and trivially on[l,l+1][l,l+1]|ψ⟩ is= 1̂ ⊗ · · · ⊗ Ô2 ⊗ · · · ⊗ 1̂. The expectation value ⟨ψ|Ô2
the other sites, Ô231

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
expressed in a simple equation by similar calculations as in the previous one:[l,l+1]=⟨ψ|Ô2
∑

|ψ⟩
(
)
(
)∗
[l] [l+1]jl+1 [l+1][l]il[l] [l+1]il+1 [l+1][l−1] [l]jl
λ[l−1]
Γ
λ
Γ
λ
⟨j
j
|
Ô
|i
i
⟩
λ
Γ
λ
Γ
λ
.
l
l+1
2
l
l+1
αl−1 αl−1 αl αl αl αl+1
αl+1
αl−1 αl−1 βl βl βl αl+1
αl+1

il ,il+1 ,jl ,jl+1
αl−1 ,αl ,αl+1
βl(3.44)
A diagrammatic proof for this formula is shown in Fig. 3.8. Again, due to the canonicalconditions of MPS, we can calculate the expectation value by using only several tensorsof MPS. Expectation values for n-site operators can be expressed in the same way.

3.3.5

Example of MPS

We show some examples of MPS representation.
S = 1/2 Néel stateS = 1/2 Néel state is the ground state of anti-ferromagnetic Ising chain and written as|Néel⟩ = |· · · ↑↓↑↓ · · ·⟩. It is a product state, and therefore has a simple MPS representation with bond dimension χ = 1:λ[2n] = λ[2n+1] = 1,
Γ[2n]↑ = 1, Γ[2n]↓ = 0,
Γ[2n+1]↑ = 0, Γ[2n+1]↓ = 1.

(3.45)(3.46)(3.47)The index αn is omitted because χn = 1 for all sites n.
S=1/2 dimer stateS = 1/2 dimer state is a chain of singlet pairs √12 (|↑⟩ |↓⟩ − |↓⟩ |↑⟩),
1
1|Dimer⟩ = √ (|↑⟩1 |↓⟩2 − |↓⟩1 |↑⟩2 ) ⊗ √ (|↑⟩3 |↓⟩4 − |↓⟩3 |↑⟩4 ) ⊗ · · · .
22(3.48)MPS for the dimer state has the bond dimension χ2n+1 = 2 and χ2n = 1:(λ[2n+1] =

√1
2

0
λ

[2n]

0
√1
2

)

(
)
(
)
, Γ[2n+1]↑ = 1 0 , Γ[2n+1]↓ = 0 1 ,

= 1,

(
Γ

[2n]↑

=

32

)
0
,
−1

[2n]↓

Γ

( )
1
=
.
0

(3.49)(3.50)

Chapter 3. Numerical Methods!
" O1[l ] " =

!

5

O.

!:.;A ":.; !:$;A":.;

":4<.; !:4;A ":4;
":4<.;

!:$;

5
5

!:9<.; ":9<.; !:9;

":4<.; !:4; ":4;!:.; ":.; !:$;

":4;

!:4;

!:9<.;A ":9<.; !:9;A
":9<.;

!:9<.;

O.
"$#":4<.; !:4; ":4;

":4<.; !:4; ":4;

"9<$#

5 5

O.

O.

"9<$#

"$#

Figure 3.7: Expectation values of one-site operators.
!
" O2[l,l +1] " =!:.; ":.;

!

":4<.; !:4; ":4; !:4C.; ":4C.;

5

":9<.; !:9;

O$

!:.;A ":.;
":.;

":4<.;!:4;A ":4; !:4C.;A":4C.;
!:4; ":4; !:4C.; ":4C.;

":4<.;

5

":9<.; !:9;A
":9<.;

O$
":4<.; !:4; ":4; !:4C.; ":4C.;

5 5

O$

Figure 3.8: Expectation values of two-site operators.

33

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
S=1/2 superposed dimer stateS = 1/2 superposed dimer state is the superposition of two dimer states,
11|Superposed Dimer⟩ ∝ · · · ⊗ √ (|↑⟩1 |↓⟩2 − |↓⟩1 |↑⟩2 ) ⊗ √ (|↑⟩3 |↓⟩4 − |↓⟩3 |↑⟩4 ) ⊗ · · ·22
1
1+ · · · ⊗ √ (|↑⟩0 |↓⟩1 − |↓⟩0 |↑⟩1 ) ⊗ √ (|↑⟩2 |↓⟩3 − |↓⟩3 |↑⟩2 ) ⊗ · · · .
22This state can be expressed as translational invariant MPS,






1/2
0√
0
0 0 0
0 1 0
λ[n] =  0 1/ 2 0  , Γ[n]↑ = 1 0 0 , Γ[n]↓ = 0 0 1 ,0 −1 0
0 0 0
0
0
1/2

(3.51)

although it is not normalized and canonicalized. This is one of the non-trivial examplesof MPS.

3.4

Time-Evolving Block Decimation method

In this section, we review the Time-Evolving Block Decimation (TEBD) method. It was
proposed by Vidal in 2004 [15] and extended to infinite systems in 2007 by himself [52].
TEBD is an algorithm to calculate imaginary or real time-evolution of a state,|ψ(τ )⟩ =e−τ H |ψ0 ⟩, |ψ(t)⟩ = e−itH |ψ0 ⟩ ,∥ e−τ H |ψ0 ⟩ ∥(3.52)

in an eﬃcient and clear way using MPS. TEBD is mathematically equivalent to well-known
DMRG method, if we formulate DMRG in MPS language. However, practically source andtendency of error in computation are diﬀerent from each other due to the diﬀerent schemesof truncation (approximation) in the algorithms. The numerical diﬀerence between TEBD
and DMRG is reviewed in Ref. [36].
There are two important points in the TEBD algorithm: decomposing time-evolution
operator into local mutually-commuting terms and maintaining canonical forms of MPS
after the action of time-evolution operator.

3.4.1

Suzuki-Trotter decomposition for time-evolution operator

We introduce Suzuki-Trotter decomposition for time-evolution operator. This decomposition makes time-evolution operator into product of local and mutually-commuting
terms.
Consider Hamiltonian that can be written as a sum of two-site operators,
H=∑nHn,n+1 =∑H2n−1,2n +

n

∑
n

34

H2n,2n+1 =: Heven + Hodd ,(3.53)

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
where Hn,n+1 is a two-site operator on site n and n + 1. First we divide time-evolution
operator9 Û (τ ) = e−τ H into small time steps, δτ = τ /N ≪ 1 (N is the total number oftime steps):()NÛ (τ ) = e−τ H = e−δτ H .
(3.54)
Then e−δτ H is decomposed into terms which contain only Heven or Hodd :e−δτ H = e−δτ (Heven +Hodd )e−δτ Heven e−δτ Hodd + O(δτ 2 ),e− δτ2 Heven e−δτ Hodd e− δτ2 Heven + O(δτ 3 ),=(1−θ)δτθδτ
e− 2 Heven e−θδτ Hodd e− 2 Heven e−(1−2θ)δτ Hodd
 · e− (1−θ)δτ
Heven −θδτ Hodd − θδτ2ee 2 Heven + O(δτ 5 ), θ = 2−211/3 ≈ 1.351.
These decompositions in each line are called 1st, 2nd and 4th order Suzuki-Trotter decomposition, respectively.
Each term of Hodd (Heven ) is mutually-commuting and the time-evolution operatore−δτ Hodd (e−δτ Heven ) can be written as a product of local terms,∏∏e−δτ Hodd =e−δτ H2n,2n+1 =:V̂2n,2n+1 .
(3.55)
n3.4.2

nUpdating canonical MPS after the action of local operator

As explained in the previous section, the canonical conditions of MPS are crucial tocalculate physical quantities eﬃciently. We need to make the canonical conditions of
MPS to hold after the time-evolution. Here we discuss the way to calculate canonical
MPS after the action of local operators. For local operators, we have to update only
several tensors to retain a canonical form of MPS.
Action of one-site unitary operator
Consider a one-site unitary operator on site l,∑Û =Uil jl |il ⟩ ⟨jl | .

(3.56)

il ,jlThe problem is how to obtain canonical MPS for |ψ ′ ⟩ = Û |ψ⟩ when
|ψ⟩ =

d
∑

χi
∑

[2]
[L−1] [L]iL
1 [i1 ] [2]i2Γ[1]iα1 λα1 Γα1 α2 λα2 · · · λαL−1 ΓαL−1 |i1 i2 · · · iL ⟩(3.57)

i1 i2 ···iL α1 ···αL−1

is canonical MPS. The answer is simple because of the locality and unitarity of the operator. We have only to change the tensor Γ[l] into∑ll(3.58)
Uil jl Γ[l]jΓ̃[l]iαl−1 αl ,αl−1 αl =jland keep the other {Γ, λ} the same. The unitarity of the matrix U ensures that the
canonical conditions hold for the updated MPS.
9When τ is real, Û (τ ) corresponds to imaginary time-evolution. For real time-evolution, we take
τ = it.

35

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

Figure 3.9: Updating MPS after the action of one-side operator Û and two-site operatorV̂ .

Action of two-site unitary operator
Consider a two-site unitary operator on site l and l + 1,
∑

V̂ =

ii

l+1
|il il+1 ⟩ ⟨jl jl+1 | .
Vjlljl+1(3.59)

il ,il+1 ,jl ,jl+1

Canonical MPS for |ψ ′ ⟩ = V̂ |ψ⟩ is obtained by modifying the tensors Γ[l] , Γ[l+1] and λ[l] .
To do this, we define a rank-2 tensor Θ(il αl−1 ),(il+1 αl+1 ) and apply SVD to it:∑ (Θ(il αl−1 ),(il+1 αl+1 ) :=)il il+1[l]jl[l] [l+1]jl+1 [l+1]λ[l−1]
ΓλΓλαl−1 αl−1 αl αl αl αl+1
αl+1 Vjl jl+1 ,(3.60)

jl ,jl+1 ,αlSVDΘ(il αl−1 ),(il+1 αl+1 ) =

χ̃l
∑

[l]
U(il αl−1 )β λ̃β Vβ(il+1 αl+1 ) =:

β=1

χ̃l
∑

[l]i

[l] [l+1]i

l+1 [l+1]l
λ[l−1]
αl−1 Γ̃αl−1 β λ̃β Γ̃βαl+1 λαl+1 .

β=1(3.61)
Note that the updated tensors can have larger bond dimension χ̃l than before. If χl−1 =χl = χl+1 = χ, the size of tensor Θ(il αl−1 ),(il+1 αl+1 ) is dχ × dχ, so χ̃l could be dχ at most.
The growth of the bond dimension is physically interpreted as increase of entanglement
between sites [1 · · · l] and sites [l + 1 · · · L] due to the action of V̂ . This is an important
point in practical numerical calculations, as we will see below.
These procedures to update MPS are summarized in Fig. 3.9.

3.4.3

TEBD algorithm

TEBD algorithm is as follows (we take imaginary time-evolution as an example. Real
time-evolution can be done in the same way):(a) prepare canonical MPS for |ψ0 ⟩ on which time-evolution operator to act.
36

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods(b) divide the time evolution operator Û (τ ) = e−τ H into small time steps and performSuzuki-Trotter decomposition. Take 1st order decomposition for example, e−δτ H =e−δτ (Heven +Hodd ) ≈ e−δτ Heven e−δτ Hodd .
(c) calculate |ψ ′ ⟩ = e−δτ Hodd |ψ0 ⟩ in MPS language. Now the time-evolution operator
is decomposed into a product of local operators V̂2n,2n+1 (Eqn. (3.55)), we should
update the tensors Γ[2n] , Γ[2n+1] and λ[2n] for every n according to Eqn. (3.61) (notethat the computations can be done in parallel for n). However, two problems appear
in the updating processes: unitarity of the time-evolution operator and truncationof the state. We discuss these two problems below.
(d) calculate e−δτ Heven |ψ ′ ⟩ in the same way and repeat the time-evolution until desired
result is obtained. In the case of imaginary time-evolution for ground state search,
convergence of energy (expectation value of Hamiltonian) in the time-evolution is
one of the indicators to stop the evolution.
Unitarity of imaginary time-evolution
Imaginary time-evolution operator Û (τ ) = e−τ H and the decomposed operator V̂2n,2n+1 =e−δτ H2n,2n+1 are not unitary when τ ̸= 0. As a result, the updated tensors Γ̃[2n] , Γ̃[2n+1] , λ̃[2n]in Eqn. (3.61) do not∑satisfy canonical conditions of MPS. Furthermore, the normalization
condition of a state χα2n
(λα2n )2 = 1 is violated after the imaginary time-evolution. In
2n =1
the TEBD algorithm for imaginary time-evolution, we have to recover the normalizationof a state by modifying Schmidt values1λ̃[2n]
λ̃[2n]→vα ,αu χ̃2n (u∑ [2n] )2t
λ̃α(3.62)

α=1

after SVD. On the other hand, the problem of canonical conditions is ignored10 in thealgorithm because the violation is small enough when δτ ≪ 1.
Truncation
In practice, we have a limit χmax for the bond dimension of MPS because an amount ofmemory storage of computers is finite. However, as we discussed, the bond dimensionof the updated tensors could be larger than before. If new bond dimension χ̃l becomeslarger than χmax , we have to truncate Schmidt vectors which have small values:′|ψ ⟩ =

χ̃l
∑

αl =1

αl =1
10

∑

χmax
[1···l]⟩≈⟩ |τ̃α[l+1···L]λ̃[l]αl |τ̃αll

Detailed discussion on this point is found in Ref. [53].

37[1···l]⟩.
⟩ |τ̃α[l+1···L]λ̃[l]αl |τ̃αll(3.63)

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
In addition, we should modify λ̃[l]1
λ̃[l] (α = 1, · · · , χmax )
λ̃[l]
α → v
uχmax ( ) αu ∑ [l] 2t
λ̃α(3.64)

α=1

in order to normalize |ψ ′ ⟩. Truncation error ϵ is defined as
(

∑dχ
ϵ=1−

α=χmax +1

[2n]

λ̃α
)2

∑dχ ( [2n]α=1 λ̃α)2.

(3.65)The truncation error gives a good estimate for accuracy of numerical calculations.
Source of error
There are two sources of error in TEBD: one from Suzuki-Trotter decomposition andthe other from truncation. The first one grows linearly in time. It can be reducedby taking smaller time step δτ or employing higher-order Suzuki-Trotter decomposition.
Both treatments increase the total number of SVD to simulate the same duration oftime. The error from truncation grows exponentially in time, and become important inlonger time scales. It is reduced by using lager χmax , which results in the increase ofcomputational cost of SVD at each time step (computational cost of SVD for dχ × dχ
matrix is O ((dχ)3 )). The characteristics of these two errors are thoroughly discussed inRef. [54].

3.5

iTEBD

In this section, infinite Time-Evolving Block Decimation (iTEBD) method is reviewed.
This method is an extension of TEBD to infinite systems. The iTEBD can simulate
systems in the thermodynamic limit directly, without any extrapolation of system size.
In this study, we use iTEBD to simulate the dynamics after the flux quench.

3.5.1

infinite MPS

In iTEBD, we assume spatial periodicity of tensors in MPS, and quantum states are
expressed by such MPS of infinite length, infinite MPS (iMPS). By assuming periodicity,
we have only to treat several tensors to simulate infinite systems. If the support of eachterm in Hamiltonian is m sites at most, the periodicity must be set to be larger thanm in order to construct iTEBD method. We note that assuming periodicity to states
is potentially dangerous because it restricts possible orders to be commensurate to that
periodicity. For example, dimerized order (which has two-site translational symmetry)cannot be expressed by iMPS with odd-number site periodicity. We should be careful
about what orders could appear when assuming the periodicity of tensors in MPS.
38

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods
The canonical conditions of iMPS are defined as the same as MPS. In the same manner
as finite MPS, the canonical conditions allow us to compute physical quantities per siteby using only a few tensors in iMPS.

3.5.2

Algorithm

The algorithm of iTEBD is simper than that of TEBD. For simplicity, consider iMPS
with two-site periodicity and Hamiltonian which consists of nearest-neighbor interaction
terms. We follow the notation of Vidal [52] below. Tensors of iMPS are written asΓ[2n] = ΓA , λ[2n] = λA , Γ[2n+1] = ΓB , λ[2n+1] = λB ,and a state |ψ⟩ is expressed as
∑
B Ai2
A Bi3
B
1
|ψ⟩ =
· · · ΓBi
α0 α1 λα1 Γα1 α2 λα2 Γα2 α3 λα3 · · · |· · · i1 i2 i3 · · ·⟩ .

(3.66)(3.67){i}{α}

Hamiltonian h of the system is∑∑∑h=h[n,n+1] =h[2n−1,2n] +h[2n,2n+1] ,nn(3.68)n

and we perform 1st order Suzuki-Trotter decomposition for the time-evolution operator
U (δτ ) := e−δτ h ,[n,n+1]:= exp(−δτ h[n.n+1] ),⊗ [2n,2n+1]⊗ [2n+1,2n+2]
AB
BA
Uδτ
:=
Uδτ
, Uδτ
:=
Uδτ
,
Uδτ

n

(3.69)(3.70)n

BA AB
U (δτ ) = UδτUδτ + O(δτ 2 ).

(3.71)BA
AB BA
3
For 2nd order decomposition, we can use U (δτ ) = Uδτ/2 Uδτ Uδτ /2 + O(δτ ).
At each time step, we update tensors ΓA , ΓB , λA (ΓA , ΓB , λB ) after the action of thedecomposed time-evolution operator U AB (U BA ). For example, after the action of U AB ,we obtain Γ̃A,B and λ̃[A] by)∑ (ij
B Ai′ A Bj ′ B
Θ(iα),(jγ) :=λα Γαβ ′ λβ ′ Γβ ′ γ λγ UiAB
(3.72)′j′i′ ,j ′ ,β ′χ̃∑Θ(iα),(jγ) =

χ̃
∑
A

A

SVD

X(iα)β λ̃A
β Yβ(jγ) =:Ai A Bj BλBα Γ̃αβ λ̃β Γ̃βγ λγ .

β=1

β=1

A diagrammatic expression of the update is shown in Fig. 3.10.

39(3.73)

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 3. Numerical Methods

Figure 3.10: Diagrammatic representations of the iTEBD algorithm. Left panel shows
the whole picture of update of tensors after the action by U AB . Right panel shows thatthe update procedure to obtain new tensors ΓA , ΓB , λA . Figures are taken from Ref. [52].

40《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4
Results and Discussions
In this chapter, we present numerical results of the dynamics of spin current and discuss
them.

4.1

Numerical data of the dynamics

In numerical calculations, we first prepare the ground state of the XXZ chain with flux)
∑ (11 −iθ − +iθ + −z z
e Si Si+1 + e Si Si+1 + ∆Si Si+1
(4.1)
Hθ = −
2
2i
by imaginary time-evolution using iTEBD method:e−τ Hθ |ψini ⟩,τ →∞ ∥ e−τ Hθ |ψini ⟩ ∥|θ⟩GS = lim

where |ψini ⟩ is initial state of the imaginary time-evolution.
Then, real time-evolution driven by the XXZ Hamiltonian without flux
)
∑ (1
1 − +
+ −
z z
H0 = −S S + S S + ∆Si Si+1
2 i i+1 2 i i+1
i

(4.2)

(4.3)

is performed on the obtained ground state |θ⟩GS ,|ψ(t)⟩ = e−iH0 t |θ⟩GS .

(4.4)

We use iTEBD again for the real time-evolution.
Numerical calculations were done in the following parameters:∆ = ±0.3, ±0.5, ±0.8, −1.0, −1.2, −1.5, −2.0, (±0.1),
π
π
π
π π π
θ = − ,− ,− ,− ,− ,− .
2
3
6 10 20 30∆ is the strength of interactions and |θ| represents the size of quench (for larger |θ| theinitial state becomes more excited one). Since we employ iTEBD, we can treat the systemin the thermodynamic limit L → ∞ directly without resorting to any extrapolation. There41

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions
is no finite-size eﬀect in our calculations. Instead the bond dimension χ of MPS determines
accuracy of calculations. As we discussed in the last chapter, in practice the truncation
error could appear in the iTEBD algorithm due to the finiteness of χ. In order to reduce
the truncation error and simulate the dynamics for longer time, we have to take larger χ.
We carefully check the truncation error and other potential errors in iTEBD calculations,and adopt data in reliable time scale. Technical details are explained in appendix A.
Numerical results of the dynamics
We show numerical results of the dynamics of the spin current J(t) in Fig. 4.1 and 4.2.
We summarize the dynamics as following. For large |θ|, θ = −π/2 and −π/3, the spincurrent shows oscillation and decay for all ∆ (the data for ∆ = ±0.1 have too long
periods of oscillations and we cannot see them clearly). For smaller |θ|, θ = −π/6 to−π/30, we observe qualitatively diﬀerent dynamics for ∆ > 0 and ∆ < 0. For ∆ > 0, theoscillation in long time scale (t ≳ 10) is not visible and the decay of the current is small(we will discuss the visible oscillation in short time scale (t ≲ 10) later). The spin current
becomes stationary within the accessible time scale of numerical calculations. For ∆ < 0,
the oscillation in long time scale is visible, although for small |θ| and |∆| the period is toolong for us to confirm the oscillation. The decay of the spin current is larger than that of∆ > 0 cases. Besides, the relaxation time is larger than in ∆ > 0 cases, and we cannot
observe stationary states withn the accessible time scale of numerical calculations.
Regarding the long-time limit of the spin current, we observe that it is non-zero ingapless case (−1 ≤ ∆ < 1) and zero in gapped case (∆ < −1) for any θ. However, we
mention that for ∆ = −0.8 and −1.0 it is diﬃcult to judge whether the long-time limit
is zero or non-zero from the numerical data because of the long time scale of relaxation.
We will discuss these points in section 4.2.
Before analyzing the numerical data quantitatively, we comment on the oscillation in
short time scale (t ≲ 10) that is visible even in ∆ > 0 for small |θ| (Fig. 4.2). This oscillation survives only immediately after the quench and it is observed for almost all parameterregions. Therefore we argue that it dose not contain much meaningful information of thesystem. We do not discuss this oscillation further in this study.
Fitting for the numerical data. To quantify characteristics of the dynamics fromnumerical data, we exploit a simple and empirical fitting function,f (t) = c + (A + B cos(ωt + ϕ)) e−t/τ , g(t) = c + Ae−t/τ .

(4.5)

We use f (t) for fitting when the oscillation in long time scale is visible. Otherwise weuse g(t) instead. Therefore, by this fitting we analyze the oscillation in long time scale
only and the oscillation in short time scale is not treated. The fitting is well for shorttimes, roughly up to the first cycle of oscillation in the case of f (t). Examples of fitting
are shown in Fig. 4.3.
Among the parameters of this fitting, we focus on c, ω and τ as the long-time limit of
the current, the frequency of oscillation and the relaxation time, respectively. Below we
discuss these three features of the dynamics in order.
42

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions

0.3

0.35

θ = - π/2

0.3∆= -0.1

J(t)

0.15

-0.5
0.5

0.1

-0.8

0.2

0.1

-0.3

0.2

J(t)

0.25

0.3

-1.0

0

-0.05
0

0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
-0.02

5

10

15
t

20

25

30

0
0.1

0.3
0.5
0.8
-0.5

5

10 15 20 25 30 35 40 45
t

0.08

-0.3

0.06
-0.8

J(t)

J(t)

0.1

0

-2.0 -1.2 -1.5

-0.05

θ = - π/6

0.04

θ = - π/10

0.02
-1.0

0

-1.2
-2.0 -1.5

-0.02
0

20

40

60

80

100

t

0.05

θ = - π/20

0.04

J(t)

0.03
J(t)

0.15

0.05

0.8

0.05

θ = - π/3

0.25

0.02
0.01
0
-0.01
0

20

40

60
t

80

100

120

0 10 20 30 40 50 60 70 80 90 100
t
0.035
0.03
0.025
0.02
θ = - π/30
0.015
0.01
0.005
0
-0.005
-0.01
0
20
40
60
80 100 120
t

Figure 4.1: Numerical data of the dynamics for θ = π/2 to−π/30.

43

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions

0.16

0.1

θ = - π/10

0.3

0.098

0.155

0.096

0.5

0.15

J(t)

J(t)

0.094

θ = - π/6

0.145
0.14

0.092
0.09
0.088

0.8

0.135
-0.5

0.13
0

0.086

-0.3

20

40

60

80

100

0 10 20 30 40 50 60 70 80 90 100
t
0.034

t

0.05

θ = - π/30

0.03350.049

0.033
0.0325

θ = - π/20
J(t)

J(t)

0.048
0.047

0.032
0.0315
0.031

0.046

0.0305
0.045

0.03
0

20

40

60
t

80

100

120

0

20

40

60
t

80

100

120

Figure 4.2: Close-up of the previous figures to see the stationary states and the oscillationin short time scale.

∆= 0.5∆= 0.8∆= -0.5∆= -0.8

θ = - π/6, line: fitting, dot: Data∆= -1.0∆= -1.2∆= -2.0

J(t)

J(t)

θ = - π/2, line: fitting, dot: Data
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
-0.05
0

5

10

15

20

25

t

0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
-0.02
0

20

40

60

80

t

Figure 4.3: Fitting for θ = −π/2 and θ = −π/6. For short times, fitting is well.

44

100

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions
1-J(t=∞)/D|θ|J(t=∞)/D

0.6

1

∆=0.3
∆=0.5
∆=0.8
∆= -0.3∆= -0.5∆= -0.8∆= -1.0
LR theory

0.5
0.4
0.3

0.1

0.2

∆=0.3
∆=0.5
∆=0.8
∆= -0.3∆= -0.5∆= -0.8∆= -1.0

0.01

0.1
0
0

0.1

0.2

0.3

0.4

0.5

0.001

0.6

0.1

|θ|

1

|θ|

Figure 4.4: Left panel: J(t = ∞)/D. Right panel: 1 − J(t = ∞)/(D|θ|). Note that wecannot obtain reliable data for ∆ = −0.5, −0.8, −1 when |θ| = π/30 ≈ 0.1.

4.2

Long-time limit of the current

As we discussed in section 2.2.1, linear response (LR) theory predicts that the long-time
limit of the current J(t = ∞) is given by the Drude weight D,
J(t = ∞) = D|θ|.

(4.6)

Exact value of the Drude weight for the XXZ chain is known as [23]D=

π sin µ, µ = arccos(−∆).
4 µ(π − µ)(4.7)

D is non-zero only for −1 ≤ ∆ < 1 and otherwise zero. Furthermore, D is even functionof ∆ except1 at ∆ = ±1. Therefore LR theory predicts that J(t = ∞) is non-zero onlyfor −1 ≤ ∆ < 1 and determined by the absolute value of ∆ (irrespective of the sign).
As mentioned in the last section, for any θ we observe that J(t = ∞) remains non-zerofor gapless case (−1 ≤ ∆ < 1) and decays to zero for gapped case (∆ < −1). Thus ournumerical data are consistent with LR prediction, as far as whether the long-time limitof the current vanishes or not.
However, the value of J(t = ∞) itself deviates from LR prediction, especially forlarge |θ|. This behavior is natural because LR theory is not reliable when |θ| is large.
We plot J(t = ∞)/D in the left panel of Fig. 4.4. If LR theory is exact, J(t = ∞)/D
must be |θ|. As expected, when |θ| becomes small J(t = ∞)/D gets to close to |θ|.
Furthermore, we find that the deviation from LR theory greatly depends on the sign of∆. To see this clearly, the normalized deviation from LR theory 1 − J(t = ∞)/D|θ|is plotted in the right panel of Fig. 4.4. For ferromagnetic interactions (∆ > 0), thenormalized deviation from LR theory shows power-law in |θ|. In contrast, for large antiferromagnetic interactions (∆ < 0), J(t = ∞) still deviates from LR theory even for small|θ|, about π/20 ≈ 0.16. We expect that if |θ| is made to be smaller than the values of ournumerical data (|θ| < π/20) the deviation will obey power-law in |θ| even for ∆ < 0 case.
1

At ∆ = 1, the Drude weight is not continuous with respect to ∆. lim D(∆) = 1/4, D(∆ ≥ 1) = 0.
∆↗1

45

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions
2.5

θ = -π/2
θ = -π/3θ = -π/6
ω = |∆|/2

2

ω

1.5
1
0.5
0
-2

-1.5

-1

-0.5

0

0.5

1

∆
Figure 4.5: Frequency ω obtained from fitting of the numerical data.

It is surprising that the deviation from LR theory, or nonlinear response of the system, islarger for ∆ < 0 than for ∆ > 0.
We may attribute the diﬀerent behavior of the deviation for ∆ > 0 and ∆ < 0
to superfluidity of the system. In one dimension, superfluidity is absent in the strict
sense. However, one can observe superfluid-like response in one dimension as dynamical
properties of a system [55]. For anti-ferromagnetic interactions (∆ < 0), or repulsive
interactions in spinless fermion language, superfluidity of the system is weaker than in
the case of ferromagnetic interactions (∆ > 0, attractive in spinless fermion). Thereforefor ∆ < 0 the large portion of normal components dissipates and J(t = ∞) gets smaller.
The smaller value of the long-time limit of the current results in2 the larger deviationfrom LR theory.
Finally, we note that this deviation is not an artifact of the fitting. For ∆ < 0 andsmall |θ|, the relaxation time is so large that it is not easy to obtain reliable value ofthe long-time limit from numerical data. Nevertheless, the decay of current for ∆ < 0 is
obviously larger than that of ∆ > 0 when we look at the dynamics (Fig. 4.1). Consideringthat in LR theory J(t = ∞) depends on the absolute value of ∆ (independent on the signof ∆), we can conclude that the large deviation from LR for ∆ < 0 is not an artifact but
an intrinsic behavior of the system.

4.3

Frequency of oscillation

The frequency of oscillation is plotted in Fig. 4.5. As already mentioned in section 4.1, we
discuss the frequency of oscillation in long time scale only. We found that the frequency
ω is proportional to |∆|. Furthermore, the proportionality constant (the slope of the plot)
is diﬀerent among θ. In gapless phase, there are no time (or energy) scales other than ∆,and therefore some dependency of the frequency upon ∆ is not so surprising. However,
the explicit form ω ∝ |∆| is not trivial. Also it should be stressed that even in gapped2

Numerically we found that J(t = ∞) is always smaller than the prediction from LR theory.

46

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions
phase (∆ < −1) the same relation ω ∝ |∆| holds. In order to explain this relation, weconceive of two scenarios.

4.3.1

Equation of motion of spin current

One possible explanation for the relation ω ∝ ∆ is as below. The second derivative of theequation of motion for the spin current isd2
ˆ (t)J(t) = ⟨[iH0 , [iH0 , J]]⟩
dt2
∑
∆2
∆2
−
z
z
⟩ (t)
Si+2
Im
⟨Si−1
Si+ Si+1= − J(t) +
4
N
i
∑(
)
∆
−
+
+
−
−
+
−
z
z
z
+
Im
⟨Si+ Si+2⟩ (t) − ⟨Si−1
Si− Si+1
Si+2⟩ (t) − ⟨Si−1
Si+ Si+1
Si+2⟩ (t) + ⟨Si−1
Siz Si+1
Si+2⟩ (t) .
Ni
(4.8)
If we can neglect the last two terms in the right hand side, this equation becomes simple
one describing harmonic oscillations, and its frequency is given by ω = |∆|/2. Although
this is an attracting explanation to the relation ω ∝ |∆|, this cannot explain the diﬀerent
proportionality constants among θ; the frequency of harmonic oscillation is independent of
the amplitude (we expect that θ is related to the amplitude of oscillation). The equationof motion could help us to understand the behavior of oscillations, although its physicalmeaning and the validity for neglecting two terms are not clear.

4.3.2

Dynamics of momentum distribution nq

As another way to reveal the nature of the dynamics and oscillation, we trace timeevolution of the momentum distribution in spinless fermion picture,
1 ∑ −iqri
nq = ⟨c̃†q c̃q ⟩, c̃q = √
e
ci .
N i

(4.9)

∑
Since the spin current is given by ⟨J(t)⟩ = N −1 q sin(q)nq (t), imbalance of population
nq between q > 0 and q < 0 means non-zero current. Figure 4.6 and 4.8 show examples ofthe dynamics of momentum distribution. At first, momentum distribution is the shifted
Fermi sea with momentum-shift by |θ| (see also Fig. 2.4). After the flux is quenched tozero at t = 0, the momentum distribution starts evolving. We observed two kinds ofoscillation in momentum distribution. One of them coincides with the oscillation of thecurrent in long time scale and could explain the origin of it. Besides, qualitatively diﬀerent
dynamics between in gapless and gapped phase is observed, which gives an interpretation
to the bahavior of the long-time limit of the current.
Two types of oscillation
Two types of oscillation are as follows. One is the collective oscillation of the momentum
distribution nq in which all momenta q participate. Amplitude of the oscillation is large for47

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions
θ= -π/6, ∆= 0.5

t=0
t=3
t=5
t=10
t=15

nq

nq

θ= -π/3, ∆= -0.5
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

0.28
0.26
0.24
0.22
0.2
0.18
0.16
0.14

J(t)

0

-1 -0.8 -0.6 -0.4 -0.2

05 10 15 20 25 30

0.2 0.4 0.6 0.8
t1

Momentum q /π

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

t=0
t=7.5
t=15

0.156
0.155
0.154
0.153
0.152
0.151

J(t)

0

10 20 30 40 50
t-1 -0.8 -0.6 -0.4 -0.20

0.2 0.4 0.6 0.8

1

Momentum q /πFigure 4.6: Time evolution of momentum distribution. Left panel: θ = −π/3, ∆ = −0.5.
Right panel: θ = −π/6, ∆ = 0.5. For the left case, the dip/peak structure and its
oscillation is observed though there are no such structures and oscillations in the rightcase.
q close to the shifted Fermi point, kF′ = π/2−|θ|. This type of oscillation is observed for all
quench parameters. The other type of oscillation is the oscillation of dip (peak) structure
far from the shifted Fermi point. Almost only two momenta pθ and pθ − π participate
in this oscillation. The visibility of the second type of oscillation almost coincides withthe existence of the oscillation in long time scale (compare two panels in Fig. 4.6). In
addition, the period of the oscillation in dip (peak) structure matches with the period of
the oscillation of the current. Therefore we expect that the dip (peak) structure is related
to the oscillation in long time scale.
Taking into account that the dip and peak are separated by momentum π, we canconstruct a physical picture of the dynamics where only two characteristic modes far
from the shifted Fermi point are scattered by the Umklapp process. The oscillation of thecurrent is governed by the oscillation of population between the two modes. This pictureis highly non-trivial from the viewpoint of universal TLL, because in TLL the low-energy
excitations nearby the Fermi points dictate all of the physics in the system. The validity of
TLL picture is ensured by equilibrium renormalization group method and of course there
are no guarantees for TLL to be valid in nonequilibrium situations. Besides, we haveto treat the shifted Fermi sea in the flux quench, which is not considered by standard
TLL theory. Previous studies [19, 20] on the interaction quench in the XXZ chain (seesection 2.1) showed that the dynamics after the interaction quench is well described by
TLL picture. As far as we know, our study is the first clear illustration in that non-TLL
structure strongly aﬀects the dynamics of the system.
Let us discuss the relationship between the dip (peak) structure and the frequency ofoscillation. From the picture where only two modes govern the dynamics of the system,
we can consider a simplified two-mode Hamiltonian
)
)(
(
) εpθ
( †
β
cpθ
†
,
(4.10)
Heﬀ = cpθ cpθ −πcpθ −πβ ∗ εpθ −πwhere εk := − cos(k) is energy dispersion relation of spinless fermions, pθ is momentumof the particle selectively scattered by the Umklapp process and β is an amplitude of the48

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions
0.5

θ= -π/2, t=6
θ= -π/3, t=9
θ= -π/4, t=15
θ= -π/6, t=30
θ= -π/10, t=30

Dip momentum pθ /πnq∆= -0.5
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

-1 -0.8 -0.6 -0.4 -0.20

0.2 0.4 0.6 0.8

0.45
0.4
0.35
0.3
0.25
0.2
-1.6

1

∆= -0.3∆= -0.5∆= -0.8

-1.4

-1.2

-1

-0.8

-0.6

-0.4

-0.2

θ

Momentum q /πFigure 4.7: Left panel: comparing the momentum distribution with diﬀerent θ at ∆ =−0.5. Right panel: the momentum of dip, pθ , versus θ.
Umklapp scattering expected to be proportional to |∆|. From this simplified two-mode
model we can deduce the frequency of the current,√ω ∝ ε2pθ + |β|2 .
(4.11)
This would conclude ω ∼ |∆| for small εpθ .
Finally, we mention the behavior of pθ as varying ∆ and θ. We calculated the valueof pθ from numerical data as the momentum at the bottom of dip structure, and plot pθ
versus θ in Fig. 4.7. Numerically we foundpθ =π− a|θ|, a ≈ 0.40 ∼ 0.53 (depending on ∆).
2(4.12)We have not yet understood the origin of this characteristic momentum pθ .
Qualitatively diﬀerent dynamics between in gapless and gapped phase
There is qualitative diﬀerence in the dynamics of momentum distribution between ingapless and gapped phase, irrespective of θ. In gapless phase, the shifted Fermi sea
structure in initial state is robust against time-evolution and it survives even after thesystem becomes stationary (Fig. 4.6). The momentum distribution of the stationary
states looks like the shifted Fermi sea again with broader shoulder than the initial states.
Shifted, imbalanced momentum distribution results in the persistent current. Therefore,
the remaining shifted Fermi sea structure agrees with the non-zero long-time limit of
the current in gapless phase. On the other hand, in gapped phase, the shifted Fermi sea
structure in initial state is not robust against time-evolution and it disappears in relatively
short time scale (Fig. 4.8). In that case the whole momentum distribution moves towards
the center (q = 0) quickly and the imbalance of nq between q > 0 and q < 0 is lost, whichbrings about J(t = ∞) = 0.
This contrasted diﬀerence in the dynamics of momentum distribution between in gapless and gapped phase might be related to additional conserved quantities recently foundonly in gapless phase [56]. Additional means that they are not obtained by the usual
algebraic Bethe ansatz. Those additional conserved quantities are responsible for ballistic transport, or non-zero Drude weight, at finite temperature in gapless phase [56]. In a49

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions
θ= -π/6, ∆= -1.2

t=0
t=1
t=2
t=3
t=10

nq

nq

θ= -π/2, ∆= -1.2
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

0.3

J(t)

0.2
0.1
0
-0.1

0 1 2 3 4 5 6 7 8 9 10

-1 -0.8 -0.6 -0.4 -0.20

0.2 0.4 0.6 0.8
t1

Momentum q /π

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

t=0
t=1
t=2
t=3
t=5
t=20
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0

J(t)

0

5

10

15 20

25

t

-1 -0.8 -0.6 -0.4 -0.20

0.2 0.4 0.6 0.8

1

Momentum q /πFigure 4.8: Time evolution of momentum distribution for gapped phase. Left panel:θ = −π/2, ∆ = −1.2. Right panel: θ = −π/6, ∆ = −1.2. In both cases, the whole
momentum distributions move towards the center (q = 0) and the current decays to zero.
similar manner, the additional conserved quantities that exist only in gapless phase would
prevents the whole movement of momentum distribution, although it is not clear how theconserved quantities aﬀect the dynamics.

4.4

Relaxation time

The relaxation time τ is plotted in Fig. 4.9. When the time scale of the decay is muchlarger than the accessible time scale of our calculations, we cannot obtain reliable τ fromfitting. We omit the data in such cases.
In general, larger |∆| and larger |θ| result in smaller τ (faster decay). This is natural
because the time derivative of the current dJ(t)/dt is proportional to ∆ (Eqn. (2.18)).
When ∆ = 0 (free fermion point), the spin current is conserved and τ must be ∞. In
agreement with this fact, τ seems to diverge as |∆| → 0. From the plots, we cannot
conclude that τ obeys exponential-law or power-law in ∆. The relationship between τand ∆ might be complicated than the simple relationship between the frequency ω and ∆,ω ∝ |∆|. We have not obtained a good explanation for the data of τ yet. In the following,
we introduce an analytical result of the relaxation time in the bosonic flux quench, as apossibly related result to our numerical data.
In bosonic case, the decay rate of the current (inverse of the relaxation time) after
the flux quench is well computed by the instanton method [57]. In that study, onedimensional Bose-Hubbard model at large filling ν is mapped to O(2) quantum rotor
model and the decay of current is treated as a quantum tunneling between states with
diﬀerent winding numbers. The decay rate of the current is estimated by the action of
instanton corresponding to such tunneling. Their conclusion is Γ ∝ p2K−2 , where Γ is the
decay rate of the current, p is momentum given to initial state and K is the Luttinger
parameter of eﬀective TLL corresponding to the Bose-Hubbard model.
Motivated by this study, we plot τ versus the Luttinger parameter K of the XXZ
Hamiltonian (K = π/(2 arccos(∆)), Eqn. (2.7)) in Fig. 4.10. The points of this plot are
not in line, which means that we cannot conclude simple relation between τ and K as in50

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 4. Results and Discussions
10000

10000

θ = -π/2
θ = -π/3θ = -π/6
θ = -π/10
θ = -π/20
θ = -π/30

1000
100

1000

τ

τ

100

10

10

1

1
0.1

0.1
-2

-1.5

-1

-0.5

0

0.5

θ = -π/10
θ = -π/20
θ = -π/30

θ = -π/2
θ = -π/3θ = -π/6

0.1

1

1

-∆

∆

Figure 4.9: Relaxation time τ . Left panel shows semi-log plot and right panel shows
log-log plot.
10000

θ = -π/2
θ = -π/3θ = -π/6
θ = -π/10
θ = -π/20
θ = -π/30

τ

1000
100
10
1
0.5

1

1.5

2

2.5

K

Figure 4.10: Relaxation time τ versus the Luttinger parameter K of the XXZ chain.
1/2 < K < 1, K = 1 and 1 < K corresponds to repulsive, free and attractive regime,
respectively.

the Bose-Hubbard model.

51

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Chapter 5Conclusion, Future WorkIn this thesis, we study the flux quench problem in the S = 1/2 XXZ chain by numericalcalculations. We focus on two aspects: the long-time limit of the spin current and its
dynamics during the time-evolution.
We find that both are dependent strongly on the anisotropy (interaction) parameter ∆of the XXZ chain and the amount of flux initially inserted. The long-time limit of thecurrent matches with the predictions of linear response (LR) theory for suﬃciently small
initial flux. However, its nonlinearity (deviation from LR theory) is largely aﬀected by
the sign of interactions. Anti-ferromagnetic interactions, or repulsive interactions for
spinless fermions, suppress superfluidity in the system and result in smaller (and also
further to LR theory’s) values of the current. Regarding the dynamics, the current decays
in time after the quench in all parameter regions and its time scale depends on thestrength of interactions and the initial flux. Furthermore, in the parameter region where
the initial flux is large or the interactions are anti-ferromagnetic (repulsive), the current
oscillates in long time scale. We numerically find that the frequency of the oscillation
is proportional to |∆|. Remarkably, the dynamics of momentum distribution of spinlessfermions indicates that this oscillation of current is caused by excitations deep inside theFermi sea. This mechanism of oscillations cannot be explained by the eﬀective Luttinger
model corresponding to the XXZ chain, which is in contrast with the previous studies on
diﬀerent types of quench in the same XXZ chain.
As future work, we must reveal the nature of specific momentum pθ that appears
in momentum distributions. It will help us to understand the oscillation of the current
completely. Besides, the origin of nonlinearity of J(t = ∞) is under investigation. It
would need deeper knowledge of superfluid-like response in one-dimensional systems. The
clear relationship between the relaxation time τ and ∆ also remains to be solved.

52

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Appendix A
Technical details on numericalcalculations
A.1

Details on numerical calculations, error estimate

For numerical simulations of the dynamics, we prepare the ground state of the XXZ
Hamiltonian with flux)
∑ (11 −iθ − +iθ + −z z
Hθ = −e Si Si+1 + e Si Si+1 + ∆Si Si+1 ,
(A.1)
2
2
i
by imaginary time-evolution of iTEBD. Second-order Suzuki-Trotter decomposition isused and typically we take χ = 500-1000. Imaginary time step δτ is reduced from δτ = 0.1
to 0.01 and 0.001 after imaginary time-evolution at coarser δτ converges. At each δτ , theconvergence is checked by the diﬀerence of energy and entanglement entropy of a half chain
between two successive time steps. We set a criterion of convergence to 1e-8 for energy
and 1e-6 for entanglement entropy. After the imaginary time-evolution at δτ = 0.001
converges, we compare the energy of the obtained states with the exact ground stateenergy of the XXZ chain [11]1 ,()∫∆ sin µ ∞tanh(x)Egs = − −dx 1 −(−1 ≤ ∆ ≤ 1)
(A.2)
4
µ
tanh(πx/µ)
0
(
)
∞
∆
1 ∑
2
Egs = − − sinh ρ+(∆ < −1)(A.3)
4
2 n=1 e2nρ + 1
where µ = arccos(−∆) and ρ = arccosh(−∆). The energy of the obtained states matcheswith the exact value for 5 or 6 digit. Therefore, even if the system is in gapless phase(−1 ≤ ∆ ≤ 1) where entanglement entropy of a half chain diverges in the thermodynamic
limit (this means χ → ∞), we confirm that the converged states of iTEBD with finite bonddimension simulate the true ground states well (related discussion is found in Ref. [20]).
1Chapter 4 in his textbook, Eqn. (4.28) and (4.44). Note that the interval of integral is misprinted
in Eqn (4.44) of his book.

53

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

6

4.5

5

4

Entanglement entropy

Entanglement entropy

Appendix A. Technical details on numerical calculations432

θ = - π/2

1

3.5
3
2.5
2
1.5

θ = - π/10

1
0.5

0
0

5

10

15
t

20

25

30

0 10 20 30 40 50 60 70 80 90 100t

Figure A.1: Entanglement entropy of a half chain after the quench. Left panel showsθ = π/2, right θ = π/10. Smaller |∆| results in slower growth of entanglement entropy.
Note that we do not show some data for the clarity of the figure.
Next, using iTEBD again, we perform real time-evolution by the XXZ Hamiltonian
without flux
)
∑ (1
1 − +
+ −
z z
H0 = −S S + S S + ∆Si Si+1 ,(A.4)
2 i i+1 2 i i+1
i
for the obtained ground states. We take time step dt = 0.01 or 0.02. As the system evolves,the entanglement entropy of a half chain grows almost linearly in time2 (Fig. A.1). This
growth of entanglement entropy increases the truncation error in the iTEBD algorithm
and it severely limits the accessible time scale of our calculations.
In order to check the validity of numerical data of the real time-evolution, three criteria
are used:(a) accumulated truncation error during the dynamics is smaller than 1e-4(b) energy conservation (unitarity of real time-evolution) is fulfilled within 1e-4∑(c) total magnetization L−1 i ⟨Siz ⟩ is zero in machine precision. We take this criterion
because initial states have no magnetization and the time-evolution by H0 preserves
the magnetization.
In addition to these criteria, for some ∆ and θ we changed the time step dt and thebond dimension χ, and confirmed the exactly same numerical data. An example of χ
dependency of the dynamics is shown in Fig. A.2.
2

In finite one-dimensional critical systems, entanglement entropy of a half chain grows linearly intime after a global quench [58]. By global we mean that the change of parameter happens homogeneously
in space.

54《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Appendix A. Technical details on numerical calculations
0.35

θ = -π/2, ∆ = -0.8

0.3

J(t)

0.25

chi=200
chi=400
chi=700
chi=1000

0.2
0.15
0.1
0.05
0
0

5

10

15
t

20

25

30

Figure A.2: An example of χ dependency of the dynamics. We adopt the data of χ = 1000
up to t = 20 (indicated by arrow).

Among these criteria, usually the one from truncation error sets the most strict limit on
the accessible time scale of calculations. In practice, especially for large |θ| and ∆ = −0.5
to −1.0, we sometimes loosened this criterion to around 1e-3. Even for such cases we
confirmed the reliability of data by checking their χ dependency or the criteria (b)(c).

A.2

iTEBD algorithm exploiting a symmetry

Considering symmetries of a system in numerical algorithms greatly reduces computational costs. As an example, let us consider exact diagonalization for a L site spin-1/2
chain. The dimension of the Hilbert space is 2L . If Hamiltonian which we want to diagonalize preserves total magnetization, or has a global U (1) (abelian) symmetry, it canbe block-diagonalized according to the value of magnetization (quantum number of thatsymmetry). The largest size of the block-Hamiltonian (the dimension of the decomposed
Hilbert space) is about L CL/2 , which is much smaller than that of the whole Hilbert space(when L = 20, 220 = 1048576 and 20 C10 = 184756). Thus computational cost of diagonalization is reduced by several magnitudes. This example clearly illustrates the advantagesof exploiting symmetries in numerical algorithms.
In tensor networks, utilizing symmetries of quantum states and Hamiltonian is also
possible [59]. For numerical calculations of this study we exploit a symmetry of the
XXZ Hamiltonian, conservation of total magnetization, to reduce computational cost3 .
In this section we describe how to incorporate that symmetry of the XXZ Hamiltonian
into the iTEBD algorithm [60]. First we prove that quantum states that have definite
magnetization can be expressed as MPS explicitly encoding the value of magnetization.
Next we describe (i)TEBD algorithm which respects the conservation of magnetization in
Hamiltonian.
3

The time needed to perform simulations was reduced by 40% in our case.

55

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Appendix A. Technical details on numerical calculations

Figure A.3: MPS which encodes magnetization. All bonds have quantum numbers andthe arrows indicate the conservation of magnetization at every site.

MPS associated with magnetization
Consider the same system as we discussed in section 3.3. The Hilbert space consists of( )⊗L
L times tensor product of local Hilbert space, H = Cd ⊗ Cd ⊗ · · · ⊗ Cd = Cd
. Forsimplicity, here we consider spin-1/2 chain case (d = 2) only. The (total) magnetizationof the system is defined as∑m̂ =Ŝiz .
(A.5)
i

When given state |ψ⟩ has definite magnetization, or being eigenstate of m̂, we canconstruct MPS for |ψ⟩ which encodes the magnetization (quantum number). In order to do
this, we must take local basis {|i⟩} as eigenstates of local magnetization: {|i⟩} = {|↑⟩ , |↓⟩}.
We associate the local basis with quantum numbers,
M (|↑⟩) = 1/2, M (|↓⟩) = −1/2,(A.6)

which is equal to the magnetization of each basis vector. Specifically, when the eigenvalueof |ψ⟩ is m, we can construct MPS
|ψ⟩ =

d
∑

χi
∑

[2]
[L−1]iL−1 [L−1] [L]iL
1 [1] [2]i2Γ[1]iα1 λα1 Γα1 α2 λα2 · · · ΓαL−2 αL−1 λαL−1 ΓαL−1 |i1 i2 · · · iL ⟩ ,(A.7)

i1 i2 ···iL α1 ···αL−1 =1

where all indices αk are associated with some quantum numbers M (αk ) and the tensors
[k]ik
Γαk−1
αk can have non-zero value only when
M (|i1 ⟩) = M (α1 )(A.8)
M (αk−1 ) + M (|ik ⟩) = M (αk ) (1 < k < L) .

M (αL−1 ) + M (|iL ⟩) = m∑[1···k][1]i [1][k]ik
Furthermore, all Schmidt vectors |ταk ⟩ = {i}{α} Γα1 1 λα1 · · · Γαk−1
αk |i1 · · · ik ⟩ become
eigenstates of magnetization, and physical meaning of M (αk ) is the eigenvalue of magne[1···k]tization for |ταk ⟩. Thus the condition (A.8) can be interpreted as the conservation ofmagnetization at site k.
Diagrammatically, the condition (A.8) is represented by drawing all bonds as arrows(Fig. A.3). The fictious arrow labeled as 0 comes in from left and the fictious arrowlabeled as m comes out to right.
56

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Appendix A. Technical details on numerical calculations
Proof∑
When |ψ⟩ =i1 ···iL ci1 ···iL |i1 · · · iL ⟩ have definite magnetization m, its wave function
ci1 i2 ···iL can be non-zero only if
M (|i1 ⟩) + M (|i2 ⟩) + · · · + M (|iL ⟩) = m.
(A.9)
(
)
Therefore, 2 × 2L−1 matrix C [1] = ci1 ,(i2 ···iL ) can be block-diagonalized according toeach value m1 := M (|i1 ⟩) by permuting the basis (i2 · · · iL ) properly. In other words, bychoosing proper 2L−1 × 2L−1 permutation matrix4 Q[1]()
A1 0 0← m1 = 1/2[1] [1]
C Q =,(A.10)
0 A2 0← m1 = −1/2
where A1 and A2 are 1 × L−1 Cm+L/2−1 and 1 × L−1 Cm+L/2 matrices. SVD for C [1] is
obtained from SVDs for A1 and A2 ,
SVD

SVD

A1 = U1 S1 V1 , A2 = U2 S2 V2 ,
(
) (
) (
)
U1 0
S1 0
V1 0 0 ( [1] )−1[1] SVD
C =:
·
·
Q
.
0 U2
0 S2
0 V2 0(A.11)

Note that unitarity of Q[1] assures that this decomposition for C [1] is actually SVD. In
the same way as explained in section 3.3, we define
(
)
(
)
(
)
U1 0
S1 0
V1 0 0 ( [1] )−1
[1]
[1]i1
[1]
Γα1 =
, λ α1 =, Vα1 ,i2 ···iL =Q
0 U2 i1 ,α1
0 S2 α1 ,α1
0 V2 0α1 ,i2 ···iLci1 ,(i2 ···iL ) =

χ1
∑

[1]

1 [1]
Γ[1]i
α1 λα1 Vα1 ,i2 ···iL .

(A.12)

α1 =1∑ [1]i[1]
From construction, Schmidt vector |τα1 ⟩ = i1 Γα1 1 |i1 ⟩ becomes an eigenvector of magnetization. We can associate index α1 with quantum number
M (α1 ) := M (|i1 ⟩) (= m − M (|i2 ⟩) + · · · + M (|iL ⟩)),

(A.13)[1]which is equal to the magnetization of Schmidt vector |τα1 ⟩.
[1][1] [1]
As a next step, we perform SVD on V α1 i2 ,i3 ···iL := λα1 Vα1 (i2 ···iL ) . From the definition[1]of V [1] , we know that V α1 i2 ,i3 ···iL can be non-zero only when
M (α1 ) + M (|i2 ⟩) + · · · + M (|iL ⟩) = m.

(A.14)[1]

Again V is block-diagonalized according to m2 := M (|α1 ⟩) + M (|i2 ⟩) by permuting
indices (α1 i2 ) and (i3 · · · iL ),
B1 0 0 0← m2 = 1[2] [1] [2]
0 B2 0 0← m2 = 0 ,P V Q =
(A.15)
0 0 B3 0← m2 = −1
when elements of matrix A = (aij ) satisfies aij = Pij , P ∈ SN (symmetric group of degree n), A iscalled permutation matrix. Permutation matrix is unitary, A† A = 1.
457

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Appendix A. Technical details on numerical calculations
where P [2] and Q[2] are 2χ1 × 2χ1 and 2L−2 × 2L−2 permutation matrices. SVDs for each
[1]
SVDmatrix Bi = Ui Si Vi yield SVD for V ,

 
 

U1 0 0
S1 0 0V1 0 0 0 (())−1[1] SVD
 0 U2 0  ·  0 S2 0  ·  0 V2 0 0 Q[2] −1 .
V=: P [2]0 0 S30 0 V3 0
0 0 U3
Then Γ[2] and λ[2] are defined by
U1 0 0
S1 0 0()()−1[1] −12
 0 U2 0 
 0 S2 0 
Γ[2]i· P [2], λ[2],
α1 α2 := λα1
α2 =0 0 U3 α1 i2 ,α2
0 0 S3 α2 ,α2V1 0 0 0 ()−1[2]Vα2 ,i3 ···iL :=  0 V2 0 0 Q[2],0 0 V3 0
α2 ,i3 ···i
L

∑
χ2

[1]

V α1 i2 ,i3 ···iL =[2][2]i2 [2]λ[1]α1 Γα1 α2 λα2 Vα2 ,i3 ···iL .

(A.16)

α2 =1

∑
[12]
[1] [2]i[1]
Schmidt vector |τα2 ⟩ = i2 α1 λα1 Γα1 α22 |τα1 ⟩ |i2 ⟩ becomes an eigenstate of magnetization,[12]and new index α2 can be associated with quantum number (magnetization of |τα2 ⟩)M (α2 ) := M (α1 ) + M (|i2 ⟩).

(A.17)

Obviously, the condition (A.8) for Γ[2] is satisfied.
[k]
Repeating this procedure (block-diagonalizing V and performing SVD for each block),we can construct MPS for |ψ⟩ exploiting quantum numbers and satisfying the condition (A.8).
■

TEBD, iTEBD algorithm with conservation of magnetization
When Hamiltonian conserves quantum numbers exploited by MPS, the block structureof {Γ[k] }k in MPS (Eqn. (A.8)) is preserved during time-evolution. In that case, TEBD
and iTEBD algorithm respecting the quantum numbers can be constructed. As we will
see later, by exploiting symmetry computational cost of the algorithm is reduced by one
or two magnitude.
Consider Hamiltonian
on spin-1/2 chain which can be decomposed as the sum of two∑site operators Ĥ = l Ĥl,l+1 and each term preserves total magnetization: [m̂, Ĥl,l+1 ] = 0.
In that case, time-evolution operators used in the (i)TEBD algorithm, such as V̂ =e−δτ Ĥl,l+1 , also preserve magnetization: [m̂, V̂ ] = 0. Hence matrix representation of V̂ =∑il il+1
ij Vjl jl+1 |il il+1 ⟩ ⟨jl jl+1 | has non-zero components only whenM (|il ⟩) + M (|il+1 ⟩) = M (|jl ⟩) + M (|jl+1 ⟩).
58(A.18)

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Appendix A. Technical details on numerical calculations
As a consequence, if initial MPS satisfies the condition (A.8), the matrix Θ considered in(i)TEBD algorithm)∑ (i il+1[l−1] [l]jl[l] [l+1]jl+1 [l+1]Θ(il αl−1 ),(il+1 αl+1 ) =λαl−1 Γαl−1 αl λαl Γαl αl+1 λαl+1 Vjlljl+1(A.19)
jl ,jl+1 ,αl

has non-zero elements only whenM (|il ⟩) + M (αl−1 ) = M (|il+1 ⟩) + M (αl+1 ).

(A.20)

Then Θ has block-diagonal structure according to ml := M (|il ⟩) + M (αl−1 ),
Θ1 0 · · ·P ΘQ =  0 Θ2 · · · ,..
.. . .
.
.
.

(A.21)

where P and Q are some permutation matrices of size 2χl−1 × 2χl−1 and 2χl+1 × 2χl+1 .
SVDs for each block Θi = Ui Si Vi yield SVD for Θ,

 
 

U1 0 · · ·S1 0 · · ·V1 0 · · ·
SVD

 
 

Θ =: P −1  0 U2 · · · ·  0 S2 · · · ·  0 V2 · · · Q−1 .
(A.22)..
.. . .
..
.. . .
..
.. . .
.
.
.
.
.
.
.
.
.
New tensors Γ̃[l] , Γ̃[l+1] , λ̃[l] are obtained by

U1 0 · · ·)−1(
 −1  0 U2 · · ·[l−1]lΓ̃[l]i,P αl−1 αl = λαl−1..
.. . .
.
.
.
iα
,α


l l−1 l


V
0
·
·
·
S
0
·
·
·
1
1
(
)−1
 0 V2 · · · −1 


[l+1]
l+1
Γ̃[l+1]i=λ, λ˜[l] αl =  0 S2 · · ·Q αl αl+1
αl+1..
.. . .
..
.. . .
.
.
.
.
.
.
α ,i
α
α ,α
l l+1 l

∑

l

l

χ̃l

Θ(il αl−1 ),(il+1 αl+1 ) =[l] [l+1]il+1 [l+1][l]il
λ[l−1]
αl−1 Γ̃αl−1 αl λ̃αl Γ̃αl αl+1 λαl+1 .

(A.23)αl =1Quantum numbers associated with new indices αl (= 1, · · · , χ̃l ) are
M (αl ) := M (|il ⟩) + M (αl−1 ) (= M (|il+1 ⟩) + M (αl+1 )) ,(A.24)∑[1···l−1][1···l][l−1] [l]il⟩ |il ⟩.
which are equal to the magnetization of Schmidt vector |ταl ⟩ = il αl−1 λαl−1 Γ̃αl−1
αl |ταl−1
By definition, the updated tensors Γ̃, λ̃ satisfy the condition (A.8).
Practical procedure of the algorithm and computational cost
In order to implement the (i)TEBD algorithm exploiting a symmetry into actual computer programs, we have to prepare lists of quantum numbers {M (αk )}k as well as MPS{Γ[k] , λ[k] }k . The update of tensors Γ[l] , Γ[l+1] is performed as below (we write bond dimension as χ):59《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Appendix A. Technical details on numerical calculations(1) computing Θi by using the list {M (αk )}k and the condition (A.20). We do not haveto calculate Θ.
(2) doing SVD for each Θi . Because the size of matrix Θi is smaller than that of Θ, thecost of SVD is greatly reduced (we will see an example below).
(3) sorting all Schmidt values obtained from the SVDs for Θi (i = 1, 2, · · · ). In practice
we merge lists of the singular value from each Θi and sort it. If the number of theSchmidt values is larger than χ in total, we have to truncate Schmidt vectors withsmaller singular values (irrespective of the number i).
(4) storing the Schmidt vectors and Schmidt values into tensors Γ̃ and λ̃ properly, byusing the list {M (αk )}k .
Let us discuss how much computational cost of SVD is reduced by exploiting the symmetry(the conservation of magnetization). SVD for n × n matrix takes O(n3 ) time. Therefore
without exploiting the symmetry, cost of SVD is O((2χ)3 ) . On the other hand, when we
consider the symmetry, we have only to perform SVD for the smaller matrix, Θi . Typically
we have experienced that the largest size of matrix Θi (i = 1, 2, · · · ) is nearly a quarter ofthe size of Θ, and therefore we gain 43 = 64 times speed-up by exploiting the symmetry5 .
For example, in the case of the ground state for ∆ = −0.5, χ = 1000, the number of indices
αk (k ∈ odd sites) that have quantum number M (αk ) = [3, 2, 1, 0, −1, −2, −3, −4, −5]
is [8, 52, 147, 241, 261, 186, 83, 20, 2], respectively6 . The largest size of Θi in this case is241 + 261 = 502, much smaller than the size of Θ, 2000.
Comment on quantum numbers in iMPS and iTEBD
In iMPS and iTEBD, the physical meaning of quantum numbers for indices αk can be
subtle because magnetization of Schmidt vector can become infinite. However, this subtlety is avoided if we consider the quantum numbers as the diﬀerence of magnetization of
Schmidt vectors between initial state and current state [60].

5We have to do many SVDs for each time step, so actual speed-up of the algorithm is less than 64
times.
6for k ∈ even sites, the number of indices αk that have M (αk ) = [7, 5, 3, 1, −1, −3, −5, −7, −9] is[2, 25, 94, 198, 267, 231, 131, 45, 7], respectively.

60

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Bibliography[1] M. Lewenstein, A. Sanpera and V. Ahufinger. Ultracold Atoms in Optical Lattices:
Simulating quantum many-body systems. Oxford University Press (2012).
[2] A. Polkovnikov, K. Sengupta, A. Silva and M. Vengalattore. “Colloquium : Nonequilibrium dynamics of closed interacting quantum systems”. Rev. Mod. Phys. 83, 863–
883 (2011).
[3] J. Dziarmaga. “Dynamics of a quantum phase transition and relaxation to a steadystate”. Advances in Physics 59, 1063–1189 (2010).
[4] M. Rigol, V. Dunjko and M. Olshanii. “Thermalization and its mechanism for genericisolated quantum systems”. Nature 452, 854–858 (2008).
[5] M. Rigol, V. Dunjko, V. Yurovsky and M. Olshanii. “Relaxation in a Completely
Integrable Many-Body Quantum System: An Ab Initio Study of the Dynamics ofthe Highly Excited States of 1D Lattice Hard-Core Bosons”. Phys. Rev. Lett. 98,
050405 (2007).
[6] M. Mierzejewski, P. Prelovšek and T. Prosen. “Breakdown of the Generalized Gibbs
Ensemble for Current-Generating Quenches”. Phys. Rev. Lett. 113, 020602 (2014).
[7] M. A. Cazalilla. “Eﬀect of Suddenly Turning on Interactions in the Luttinger Model”.
Phys. Rev. Lett. 97, 156403 (2006).
[8] N. Nessi, A. Iucci and M. A. Cazalilla. “Quantum Quench and Prethermalization
Dynamics in a Two-Dimensional Fermi Gas with Long-Range Interactions”. Phys.
Rev. Lett. 113, 210402 (2014).
[9] T. Giamarchi. Quantum physics in one dimension. Number 121 in The international
series of monographs on physics. Clarendon Press (2004).
[10] A. A. Abrikosov, L. P. Gorkov and I. E. Dzyaloshinski. Methods of quantum field
theory in statistical physics. Courier Dover Publications (1975).
[11] M. Takahashi. Thermodynamics of One-Dimensional Solvable Models. Cambridge
University Press (1999). Cambridge Books Online.
[12] P. Di Francesco, P. Mathieu and D. Sénéchal. Conformal field theory. Graduate texts
in contemporary physics. Springer (1997).
61

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Bibliography[13] A. W. Sandvik. “Computational Studies of Quantum Spin Systems”. In AIP Conf.
Proc, volume 1297, 135 (2010).
[14] S. White. “Density matrix formulation for quantum renormalization groups”. Phys.
Rev. Lett. 69, 2863–2866 (1992).
[15] G. Vidal. “Eﬃcient Simulation of One-Dimensional Quantum Many-Body Systems”.
Phys. Rev. Lett. 93, 040502 (2004).
[16] E. Lieb, T. Schultz and D. Mattis. “Two soluble models of an antiferromagnetic
chain”. Annals of Physics 16, 407 – 466 (1961).
[17] P. Barmettler, M. Punk, V. Gritsev, E. Demler and E. Altman. “Relaxation of
Antiferromagnetic Order in Spin-1/2 Chains Following a Quantum Quench”. Phys.
Rev. Lett. 102, 130603 (2009).
[18] M. Fagotti, M. Collura, F. H. L. Essler and P. Calabrese. “Relaxation after quantum
quenches in the spin- 21 Heisenberg XXZ chain”. Phys. Rev. B 89, 125101 (2014).
[19] C. Karrasch, J. Rentrop, D. Schuricht and V. Meden. “Luttinger-liquid universality
in the time evolution after an interaction quench”. Phys. Rev. Lett. 109, 126406(2012).
[20] F. Pollmann, M. Haque and B. Dóra. “Linear quantum quench in the Heisenberg
XXZ chain: Time-dependent Luttinger-model description of a lattice system”. Phys.
Rev. B 87, 041109 (2013).
[21] J.-S. Caux and F. H. L. Essler. “Time Evolution of Local Observables After Quenching to an Integrable Model”. Phys. Rev. Lett. 110, 257203 (2013).
[22] W. Liu and N. Andrei. “Quench Dynamics of the Anisotropic Heisenberg Model”.
Phys. Rev. Lett. 112, 257204 (2014).
[23] B. Shastry and B. Sutherland. “Twisted boundary conditions and eﬀective mass in
Heisenberg-Ising and Hubbard rings”. Phys. Rev. Lett. 65, 243–246 (1990).
[24] B. Sutherland and B. Shastry. “Adiabatic transport properties of an exactly soluble
one-dimensional quantum many-body problem”. Phys. Rev. Lett. 65, 1833–1837(1990).
[25] V. KOREPIN and A. WU. “ADIABATIC TRANSPORT PROPERTIES AND
BERRY S PHASE IN HEISENBERG-ISING RING”. International Journal of Modern Physics B 05, 497–507 (1991).
[26] C. J. Hamer, G. R. W. Quispel and M. T. Batchelor. “Conformal anomaly andsurface energy for Potts and Ashkin-Teller quantum chains”. Journal of Physics A:
Mathematical and General 20, 5677 (1987).
[27] J. Mun, P. Medley, G. K. Campbell, L. G. Marcassa, D. E. Pritchard and W. Ketterle.
“Phase Diagram for a Bose-Einstein Condensate Moving in an Optical Lattice”. Phys.
Rev. Lett. 99, 150604 (2007).
62

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Bibliography[28] S. Peotta, C.-C. Chien and M. Di Ventra. “Phase-induced transport in atomic gases:
From superfluid to Mott insulator”. Phys. Rev. A 90, 053615 (2014).
[29] S. Hild, T. Fukuhara, P. Schauß, J. Zeiher, M. Knap, E. Demler, I. Bloch and
C. Gross. “Far-from-Equilibrium Spin Transport in Heisenberg Quantum Magnets”.
Phys. Rev. Lett. 113, 147205 (2014).
[30] J. Schachenmayer, G. Pupillo and A. J. Daley. “Time-dependent currents of onedimensional bosons in an optical lattice”. New Journal of Physics 12, 025014 (2010).
[31] R. Orús. “A practical introduction to tensor networks: Matrix product states andprojected entangled pair states”. Annals of Physics 349, 117 – 158 (2014).
[32] M. Troyer and U.-J. Wiese. “Computational Complexity and Fundamental Limitations to Fermionic Quantum Monte Carlo Simulations”. Phys. Rev. Lett. 94, 170201(2005).
[33] P. Corboz, R. Orús, B. Bauer and G. Vidal. “Simulation of strongly correlated
fermions in two spatial dimensions with fermionic projected entangled-pair states”.
Phys. Rev. B 81, 165104 (2010).
[34] F. Verstraete and J. I. Cirac. “Renormalization algorithms for Quantum-Many Body
Systems in two and higher dimensions”. arXiv:cond-mat/0407066 (2004).
[35] G. Vidal. “Class of Quantum Many-Body States That Can Be Eﬃciently Simulated”.
Phys. Rev. Lett. 101, 110501 (2008).
[36] U. Schollwöck. “The density-matrix renormalization group in the age of matrix
product states”. Annals of Physics 326, 96 – 192 (2011).
[37] X. G. Wen. “TOPOLOGICAL ORDERS IN RIGID STATES”. International Journalof Modern Physics B 04, 239–271 (1990).
[38] X. Chen, Z.-C. Gu and X.-G. Wen. “Local unitary transformation, long-range quantum entanglement, wave function renormalization, and topological order”. Phys.
Rev. B 82, 155138 (2010).
[39] X. Chen, Z.-C. Gu, Z.-X. Liu and X.-G. Wen. “Symmetry protected topologicalorders and the group cohomology of their symmetry group”. Phys. Rev. B 87,
155114 (2013).
[40] F. Pollmann, A. Turner, E. Berg and M. Oshikawa. “Entanglement spectrum of a
topological phase in one dimension”. Phys. Rev. B 81, 064439 (2010).
[41] F. Pollmann, E. Berg, A. Turner and M. Oshikawa. “Symmetry protection of topological phases in one-dimensional quantum spin systems”. Phys. Rev. B 85, 075125(2012).
[42] J. Eisert, M. Cramer and M. Plenio. “Colloquium : Area laws for the entanglement
entropy”. Rev. Mod. Phys. 82, 277–306 (2010).
63

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Bibliography[43] M. Hastings and T. Koma. “Spectral Gap and Exponential Decay of Correlations”.
Communications in Mathematical Physics 265, 781–804 (2006).
[44] M. B. Hastings. “An area law for one-dimensional quantum systems”. Journal of
Statistical Mechanics: Theory and Experiment 2007, P08024 (2007).
[45] M. Wolf, F. Verstraete, M. Hastings and J. Cirac. “Area Laws in Quantum Systems:
Mutual Information and Correlations”. Phys. Rev. Lett. 100, 070502 (2008).
[46] P. Calabrese and J. Cardy. “Entanglement entropy and quantum field theory”. Journal of Statistical Mechanics: Theory and Experiment 2004, P06002 (2004).
[47] S. Sen. “Average Entropy of a Quantum Subsystem”. Phys. Rev. Lett. 77, 1–3(1996).
[48] D. Poulin, A. Qarry, R. Somma and F. Verstraete. “Quantum Simulation of TimeDependent Hamiltonians and the Convenient Illusion of Hilbert Space”. Phys. Rev.
Lett. 106, 170501 (2011).
[49] M. Hastings. “Solving gapped Hamiltonians locally”. Phys. Rev. B 73, 085115(2006).
[50] G. Golub and C. Van Loan. Matrix Computations. Johns Hopkins studies in themathematical sciences. Johns Hopkins University Press (1996).
[51] D. Perez-Garcia, F. Verstraete, M. M. Wolf and J. I. Cirac. “Matrix Product State
Representations”. Quantum Info. Comput. 7, 401–430 (2007).
[52] G. Vidal. “Classical Simulation of Infinite-Size Quantum Lattice Systems in One
Spatial Dimension”. Phys. Rev. Lett. 98, 070201 (2007).
[53] R. Orús and G. Vidal. “Infinite time-evolving block decimation algorithm beyond
unitary evolution”. Phys. Rev. B 78, 155117 (2008).
[54] D. Gobert, C. Kollath, U. Schollwöck and G. Schütz. “Real-time dynamics in spin- 12
chains with adaptive time-dependent density matrix renormalization group”. Phys.
Rev. E 71, 036102 (2005).
[55] T. Eggel, M. Cazalilla and M. Oshikawa. “Dynamical Theory of Superfluidity in One
Dimension”. Phys. Rev. Lett. 107, 275302 (2011).
[56] T. Prosen and E. Ilievski. “Families of Quasilocal Conservation Laws and Quantum
Spin Transport”. Phys. Rev. Lett. 111, 057203 (2013).
[57] I. Danshita and A. Polkovnikov. “Quantum phase slips in one-dimensional superfluids
in a periodic potential”. Phys. Rev. A 85, 023638 (2012).
[58] P. Calabrese and J. Cardy. “Evolution of entanglement entropy in one-dimensional
systems”. Journal of Statistical Mechanics: Theory and Experiment 2005, P04010(2005).
64

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Bibliography[59] S. Singh, R. N. C. Pfeifer and G. Vidal. “Tensor network decompositions in thepresence of a global symmetry”. Phys. Rev. A 82, 050301 (2010).
[60] M. B. Hastings. “Light-cone matrix product”. Journal of Mathematical Physics 50,
095207 (2009).

65

《修士論文》

物性研究・電子版 Vol. 4, No. 3 043601 （2015年8月）

Acknowledgement
I would like to thank my supervisor Masaki Oshikawa for his comments, remarks and suggestions throughout my master’s course. Furthermore, I express my gratitude to Grégoire
Misguich for introducing me the problem studied in this thesis. Without his kind comments, helpful discussions and encouragement, I could not complete my master’s thesis.
Discussions with Ippei Danshita, Shintaro Takayoshi and Shunsuke Furukawa were
helpful for me to carry out the study. I would like to acknowledge them.
I would like to thank our group members: Yasuhiro Tada, Wenxing Nie, Miklos Lajko, James Quach, Yuta Kumano, Yohei Fuji, Emika Takata, Soichiro Mohri, Yoshiki
Fukusumi, Hiroyuki Fujita and secretary Atsuko Tsuji. My life in Kashiwa with them
was really enjoyable and academically fruitful.
Finally, I’m happy to thank my parents for their patients and constant supports inmy whole life.
The numerical calculations of this study were done mainly in Super Computer Centerof ISSP. I acknowledge Advanced Leading Graduate Course for Photon Science (ALPS) for
financial support and my sub-supervisor in ALPS, Professor Makoto Kuwata-Gonokami.

66