2006.08.03 Computational Approaches to Quantum Critical Phenomena @ ISSP

Quantum Phase Transitions in
Low Dimensional Magnets
Synge Todo (藤堂眞治)
Department of Applied Physics, University of Tokyo<wistaria@ap.t.u-tokyo.ac.jp>

Low-Dimensional Magnets• Spin-S Heisenberg antiferromagnets

H =J

"
!

!i,j"

1 + −
− +
z z
Si Sj − (Si Sj + Si Sj )2#(non-frustrated, no external field)• 3D or higher: thermal phase transitions at finite T• 2D: Néel long-range order at ground state• 1D: Haldane conjecture• S=1/2, 3/2,... : critical (i.e. quasi-LRO) GS• S=1,2,... : singlet GS with finite gap

2D Heisenberg Antiferromagnet
H=

!
i,j

!

S2i,j · S2i+1,j + α

spatial structurei,jS2i+1,j · S2i+2,j + J

!

!
i,j

Si,j · Si,j+1

ground state phase diagram

S=1

isotropic
squarelattice

H-II
10

S=1 ladder 1

2D QCP AF

J’
J’

0.1
0.08

1

!

M.Matsumoto, et al (2002)

AF

0.06

antiferro
LRO

VBS
state H-II
0.04
0.02
0

decoupled
dimers

H-I
0

0.2

0.6
0.8
1
bond 0.4
alternating !
1D QCP
chain

Haldanechain

Agenda• (Very brief) review on loop cluster quantum

Monte Carlo• Gap estimation by using loop cluster QMC• extended moment method• Topological order in 1D magnets• twist order parameter and its improved

estimator• Quantum surface transition in 2D• perturbation by using QMC

Path-Integral Representation• H = H0 (diagonal part) + V (off-diagonal part)Z = tr e−βH0

U (β) with e−β(H0 +V ) ≡ e−βH0 U (β)• Time-dependent perturbation expansion

V (τ ) ≡ eτ H0 V e−τ H0
dU (τ )/dτ = −V (τ )U (τ )! β
U (β) = 1 −V (τ )U (τ ) dτ
!0 β
! β ! τ1
=1−
V (τ1 ) dτ1 +V (τ1 )V (τ2 ) dτ1 dτ2 + · · ·0!

Z = tr e−βH0 1 +

0

0

∞
"

n=1

(−1)n

# β # τ100···

# τn−1 $n0i=1V (τi )

n
$

i=1

dτi

%

World-Line Configuration• Each term is specified by (α, n, {bi , τi })“world-line configuration”

τ=β

τ8

τ7
b=1

b=4

τ6

τ5

τ4
b=2

τ2
b=1

b=4

τ3

n=8

|α! = | ↑↓↓↑↓↑ !

b=3

τ1
b=3

τ=0

• Weight of each configuration can be calculated easilyw(α, n, {bi , τi }) ≡ (−1)n #α|e−βH0 Vb1 (τ1 )Vb2 (τ2 ) · · · Vbn (τn )|α$Path Integral Representation• A d-dimensional quantum system is mapped

to a (d+1)-dimensional classical system.
additional dimension = imaginary time
• Up (down) spins form world lines due to Sz

conservation.
• Weight of mapped configuration (world lines)

is expressed as a product of local weights.
• Local updates or cluster updates (loopupdate) can be performed on world lines.

High Temperature Series
Representation (SSE)• Break up the Hamiltonian into diagonal andoff-diagonal bond terms−H =diagonal
off-diagonal

! !

k=1,2

Hk,b

b

h z
z z
H1,b = −Jsi sj + (si + szj )

z
J + −
− +
H2,b = (si sj + si sj )2

• Taylor expansion of density matrix (ip={k,b})

Z=

!
α

!α|e−βH |α" =

∞ ! n
!!
β
α n=0 Snn!

!α|Hi1 Hi2 · · · Hin |α"

Sandvik and Kurkijarvi (1991)

World-Line Configuration in HTSR
• Each term is specified by (α, {Hip })“operator string”

Z=

!
α

!α|e−βH |α" =

∞ ! n
!!
β
α n=0 Snn!

!α|Hi1 Hi2 · · · Hin |α"High Temperature Series
Representation (SSE)• High temperature series always converges

due to the finiteness of the system(cf. HT series in thermodynamic limit)• Similarity to path-integral representation• Integer index instead of continuous variable

in imaginary time direction

Non-Local Updates• Non-local updates are necessary for simulating in

grand canonical ensemble (different winding numberof world lines)• Strong restriction in world-line configuration:

All possible changes of world lines form “loops” dueto local Sz conservation. (Finding allowed local
update is also very tough task!)• Non-local updates may reduce auto-correlation time

loop cluster update

One MCS of Loop AlgorithmPI

breakuploop flip

HTS

diagonal update

loop flip

Loop Cluster Update• Changes shape of world lines as well as their winding

number• Does not change expansion order n in HTSR

Not ergodic by itself (used with diagonal update)• Loop cluster update for spin-S XYZ models• No loop solution in the presence of (longitudinal)

external field (need to introduce global weight to loops)• Both in PI and HTS representations

C++ implementation for S=1/2 Heisenberg antiferromagnet

path integral

high temperature seriesfor (unsigned int mcs = 0; mcs < therm + sweeps; ++mcs) {std::swap(operators, operators_p); operators.clear();fragments.resize(nsites);
std::fill(fragments.begin(), fragments.end(), fragment_t());
std::copy(boost::counting_iterator<int>(0),boost::counting_iterator<int>(nsites), current.begin());
double t = r_time();for (std::vector<local_operator_t>::iterator opi = operators_p.begin();t < 1 || opi != operators_p.end();) {if (opi == operators_p.end() || t < opi->time) {unsigned int b = nbonds * random();
if (spins[left(nbonds, b)] != spins[right(nbonds, b)]) {operators.push_back(local_operator_t(b, t)); t += r_time();} else {t += r_time(); continue;}} else {if (opi->type == diagonal) {++opi; continue;} else {operators.push_back(*opi); ++opi;}}
std::vector<local_operator_t>::iterator oi = operators.end() - 1;
unsigned int s0 = left(nbonds, oi->bond);
unsigned int s1 = right(nbonds, oi->bond);oi->lower_loop = unify(fragments, current[s0], current[s1]);
oi->upper_loop = current[s0] = current[s1] = add(fragments);
if (oi->type == offdiagonal) { spins[s0] ^= 1; spins[s1] ^= 1; }}}}for (int s = 0; s < nsites; ++s) unify(fragments, s, current[s]);int nc = 0;
for (std::vector<fragment_t>::iterator ci = fragments.begin();ci != fragments.end(); ++ci) if (ci->is_root()) ci->id = nc++;clusters.resize(nc);for (std::vector<fragment_t>::iterator ci = fragments.begin();ci != fragments.end(); ++ci) ci->id = cluster_id(fragments, *ci);for (std::vector<cluster_t>::iterator pi = clusters.begin();pi != clusters.end(); ++pi) *pi = cluster_t(random() < 0.5);
for (std::vector<local_operator_t>::iterator oi = operators.begin();oi != operators.end(); ++oi)
if (clusters[fragments[oi->lower_loop].id].to_flip ^clusters[fragments[oi->upper_loop].id].to_flip) oi->flip();for (unsigned int s = 0; s < nsites; ++s)if (clusters[fragments[s].id].to_flip) spins[s] ^= 1;}

for (unsigned int mcs = 0; mcs < therm + sweeps; ++mcs) {int nop = operators.size();
std::swap(operators, operators_p); operators.clear();fragments.resize(nsites);
std::fill(fragments.begin(), fragments.end(), fragment_t());
std::copy(boost::counting_iterator<int>(0),boost::counting_iterator<int>(nsites), current.begin());
bool try_gap = true;for (std::vector<local_operator_t>::iterator opi = operators_p.begin();
try_gap || opi != operators_p.end();) {if (try_gap) {unsigned int b = nbonds * random();
if (spins[left(nbonds, b)] != spins[right(nbonds, b)] &&(nop + 1) * random() < lb2) {operators.push_back(local_operator_t(b)); ++nop;} else {try_gap = false; continue;}} else {if (opi->type == diagonal && lb2 * random() < nop) {--nop; ++opi; continue;} else {operators.push_back(*opi); ++opi; try_gap = true;}}
std::vector<local_operator_t>::iterator oi = operators.end() - 1;
unsigned int s0 = left(nbonds, oi->bond);
unsigned int s1 = right(nbonds, oi->bond);oi->lower_loop = unify(fragments, current[s0], current[s1]);
oi->upper_loop = current[s0] = current[s1] = add(fragments);
if (oi->type == offdiagonal) { spins[s0] ^= 1; spins[s1] ^= 1; }}for (int s = 0; s < nsites; ++s) unify(fragments, s, current[s]);int nc = 0;
for (std::vector<fragment_t>::iterator ci = fragments.begin();ci != fragments.end(); ++ci) if (ci->is_root()) ci->id = nc++;clusters.resize(nc);for (std::vector<fragment_t>::iterator ci = fragments.begin();ci != fragments.end(); ++ci) ci->id = cluster_id(fragments, *ci);for (std::vector<cluster_t>::iterator pi = clusters.begin();pi != clusters.end(); ++pi) *pi = cluster_t(random() < 0.5);
for (std::vector<local_operator_t>::iterator oi = operators.begin();oi != operators.end(); ++oi)
if (clusters[fragments[oi->lower_loop].id].to_flip ^clusters[fragments[oi->upper_loop].id].to_flip) oi->flip();for (unsigned int s = 0; s < nsites; ++s)if (clusters[fragments[s].id].to_flip) spins[s] ^= 1;}

Improved Estimator• Fortuin-Kasteleyn RepresentaionZ=!

W (c) =c!

W (c, g) =c,g!

∆(c, g)f (g)c,g• In loop algorithm, graphs are generated accordingtoW (g) ≡!

W (c, g) = 2Nc f (g)c• Expected value of a quantity: A(c)!A(c)" = Z −1!

A(c)W (c)c!
!
= Z −1A(c)∆(c, g)f (g)cg!
!
= Z −1
A(g)W (g)A(g) ≡
A(c)∆(c, g)gcImproved Estimator−1!A(c)" ≈ NMC!

• e.g.) spin correlation function

A(g)Ci,j (c) = Siz Sjz• spin correlation represented in “graph language”|Ci,j (g)| =!

1/4 if i and j are on the same loop0otherwise• improved statistics (i.e. small statistical errors)• “loops” (or “graphs”) themselves have a physicalmeaning

Gap Estimation by using loop cluster QMC

How to Calculate Gap by QMC• “gap” is easily calculated in exact diagonalization∆ = E1 − E0

but not in world-line QMC• fitting temperature dependence of susceptibility, etc
0.3χ ∼ exp(−∆/T )0.25• temperature range?

0.2

~
!

• error estimation?

spin-1 ladder

0.15

0.3
0.25
0.2
0.15
0.1
0.05
0
0.01

0.1
0.05
0
0

0.5

1

~
T

0.1

1
1.5

2

A Tiny Gap in a Huge Gap• frequency distribution at∆1

T ! ∆1d0
P0 =d0 + d1 exp(−∆0 /T )• ex) edge state in open S=1 chain
.....

∆0
∆0 ! ∆1effective interaction= singlet-triplet gap << bulk gap ~ 0.41J

Imaginary-Time Correlation Function• We want a generic and systematic way of measuring

gap• Correlation function in imaginary-time direction

C(τ ) =
=

!

i,j
!

!(−1)

|i−j|

z
z
Si (0)Sj (τ )#ai cosh[(τ − β/2)/ξi ]ξ0 = ∆−1> ξ1 > · · ·• However, numerical inverse Laplace transformation is

ill-posed problem (c.f. MaxEnt)• Systematic improvement of moment method

Moment Method• Fourier transform of correlation functionC̃(ω) =! βC(τ )eiωτ dτ

0• Effectively measured by using improved estimator

τj+1
τj

#2 %
! " ##$#iωτ# Si (τ )e dτ #
=
#
#
! #
#2 $
!"
"
#
#
iω(τj+1 −τj )=Si (τj )e
/ω #
#
!

j

• Second-moment method (0-th order approximation)

β
(0)
ˆ
ξ =
2π

!

"
# ai ξi $C̃(0)− 1 → ξ0 1 −2a0 ξ0C̃(2π/β)i=1Higher-Order Approximations• Second-order estimator

β
(2)
ˆ
ξ =
4π

!

"

# ai
C̃(0) − C̃(2π/β)− 1 → ξ0 1 −32a0C̃(2π/β) − C̃(4π/β)
i=1• Fourth-order estimator

β
(4)
ˆ
ξ =
6π

$

ξi
ξ0

%3 &

!
!

3C̃(0) − 4C̃(2π/β) + 4C̃(4π/β)10−15C̃(2π/β) − 8C̃(4π/β) + 3C̃(6π/β)→ ξ0 1 −

" ai
i=1

2a0

• convergence test

#

xii
ξ0

$5 %

S=1, L=128, T=1/64(0)ˆ0.41127(2)
1/ξ
(2)
ˆ
0.41048(6)1/ξ1/ξˆ(4) 0.4104(1)ED

0.41048(2)Haldane Gap of S=3 Chain• Estimation of extremely large correlation lengthand small gap
1000

ξ = 637(1)∆ = 0.01002(3)J

Simulation

L = 5792
T = 0.001

correlation length (!x, !")S=3 Chain Results

S=1
2
3
!x
!"

100

10

1

TL=2S (fixed)
0.001

0.01

0.1

1

T

Todo and Kato (2001)

Haldane gap of S=1,2,3 chains∆

method
S=1 0.413(7)MCPM (Nightingale-Blöte 1986)
0.4150(2) 6.03(2)DMRG (White-Huse 1992)
0.41049(2) 6.2
ED (Golinelli et al 1994)
0.408(12)QMC (Yamamoto 1995)
0.41048(6) 6.0164(2) QMC+loop (Todo-Kato 2001)
S=2 0.074(16)QMC (Yamamoto 1995)
0.055(15)DMRG (Nishiyama et al 1995)
0.085(5)49(1)
DMRG (Schollwöck-Jolicœur 1995)
0.090(5)50(1)
QMC+loop (Kim et al 1997)
0.0876(13)DMRG (Wang et al 1999)
0.08916(5) 49.49(1) QMC+loop (Todo-Kato 2001)
S=3 0.01002(3) 637(1)QMC+loop (Todo-Kato 2001)

ξxFull Dispersion Curve• (In principle) moment2
1.8

S = 1/2

1.6

L = 256

1.4
!(k)

1.2
1
0.8
exact
T=0.2J
T=0.1J
T=0.05J
T=0.02J

0.6
0.4
0.2
0
0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1method can be extended
to any k-values• But (in practice) suffersfrom large systematic andstatistical errors

k/"

2

3.5
S=1

S = 1/2
L = 256

1.5

2.5
!(k)

!(k)

3

1
exact
T=0.2J
T=0.1J
T=0.05J
T=0.02J

0.5

0
0

2
1.5
1

L=16
L=32
L=64

0.5

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
k/"

1

0
0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
k/"1

Topological Order in 1D Magnets

Topological order in 1D magnets• S=1 antiferromagnetic Heisenberg ladder
0.3

J

0.25

K
0.2

~
!

J

0.15

• no QPT between two limits

0.3
0.25
0.2
0.15
0.1
0.05
0
0.01

0.1
0.05

Haldane state

0

0

0.5

1

0.1

1
1.5

~
T

2

1

dimer state

↕

~
!
0.1

"

0.01

ST, M.Matsumoto, C.Yasuda, and H.Takayama (2001)

0

0.2

#1

0.4

R=K/J
0.6

R

0.8

1

Lieb-Schultz-Mattis Theorem•

Lieb-Schultz-Mattis overlap functionzL ≡ "Ψ0 | exp[i••|Ψ0 ! : ground state
exp[i

L
!
2π

L j=1

L
!
2π

L j=1

z
jSj ]|Ψ0 #jSjz ]|Ψ0 ! : “twisted” state with

E ! E0 + O(L−1 )zL → 0 (L → ∞)gapless or degenerating ground state

E.Lieb, T. Schultz, D.Mattis 1961

Twist Order Parameter• QPT in spin-S bond-alternating chain as arearrangement of VBS pattern

S=3/2

S=1

S=1/2

• Expected value of LSM twist-

1

operator= twist order parameter
2π
zL ≡ "Ψ0 | exp[ijSjz ]|Ψ0 #L j=1
M.Nakamura and ST (2002)

0.5

(1)

S= 2
zL

L
!

S= 1

0

-0.5

-1

3
S= 2

-1

1
S= 2

-0.5

0
δ

0.5

1

Sign of Twist Order Parameter• contribution from a dimeriπcos[ |i − j|]Lj• S=1 dimer phase: (2,0) VBS state
.....
1

2

3

4

L-1

L

π L
[cos ] → +1
L

• Haldane phase: (1,1) VBS state
.....
1

2

3

zL → (−1)b

4

L-1

L

π L−1
π
[cos ][cos (L − 1)] → −1
LL

Effects of Randomness on zL• weak bond randomness

rearrangement of local VBS pattern

local rearrangement of VBS pattern can notchange the parity at boundaryzL could be a good topological orderparameter also for random systems
T.Arakawa, ST, and H.Takayama (2005)

Twist Order Parameter in
Loop Algorithm• Split-spin representation of S=1 chain
k=1
2
j=1

. . . .

2

L

• Twist order parameter

zL =

!

exp

L #
" 2πi #Lj ŝzj,k

j=1 k=1,2$%

MCS• Decomposing into contributions from each loop

exp

L "
! 2πi "

L

j=1 k=1,2#

j ŝzj,k = exp

Nc! πi "

L

!=1

#

X! =

Nc
$

exp

!=1

X! ≡

!

(j,k)∈C!

! πi
L

X!

2 j ŝzj,k#Improved Estimator• Averaging over

1 !

2 Nc

···

d1 =±1

=

Nc
!
1"

!=1

2

2Nc configurations to obtain “A(g)”

Nc
! "

exp

dN =±1 !=1

exp

# πi
L

X!

$

# πi
L

N
#π
$
πi $% !c+ exp − X! =
|X! |
cos
L
L

• Contribution from a loopcos
• for 2-site loop !

d! X!

!π

L

#

$

"
|X! |

!=1

"
"
!π
π
|X! | = cos|i − j|cos
LL= contribution from a dimer

Loops vs VBS picture• loop defines a singlet pair (or group) of spins in the

VBS picture• loop across the periodic boundary 4n (4n+2) times

contributes +1 (-1) to twist-order parameter• local fluctuation of loop structure can not change thesign (as seen in the VBS picture for random system):well-defined (topologically protected) sign in each
VBS phase.
• zL = 0 means the existence of global loop(s), i.e.

correlation length ~ system size

Spin-1 ladder• Plaquette singlet solid (PSS) state in S=1 ladder
J
K
J

Haldane state

uniform susceptibility
0.3

↕

0.2

~
!

PSS state

0.25

dimer state

↕

0.15

0.3
0.25
0.2
0.15
0.1
0.05
0
0.01

0.1
0.05
0
0

ST, M.Matsumoto, C.Yasuda, and H.Takayama (2001)

0.5

1

~
T

0.1

1
1.5

2

Quantum Surface Transition

Edge State of Open S=1 chain• Open S=1 chain with odd L
.....

1

0.7

Sz=0 and ±1

0.6

0.6

frequency

frequency

0.8
Sz=±1

0.4
0.2
0
0.0001

0.4

0.2

L = 7,9,...,33
0.001

0.01
T

0.1

Sz=0

!p = 6.0153 (Todo and Kato 2001)
0.1

1

scaling plot

0.5

0.3

Sz=0

Sz=±1

1
T exp[L/!p]

10

Jeff ≈ 0.754 exp(−L/6.01)

0-dim surface of 1-dim systemgeneralization to higher-dimensions?

Spin-1 Strips with Finite Width
H=

!
i,j

S2i,j · S2i+1,j + α!
i,jS2i+1,j · S2i+2,j + J

!

!
i,j

Si,j · Si,j+1

H-II
10

L

J’

1

1

!

J’

∞
Periodic boundary conditions
in y-direction and L → ∞

AF

0.1
0.08

AF

0.06

H-II

0.04
0.02
0

H-I
0

What happens for open boundaries?

0.2

0.4

0.6

0.8

1

!

M.Matsumoto, et al (2002)

Open S=1 Strips
H=

!
i,j

S2i,j · S2i+1,j + α!
i,jS2i+1,j · S2i+2,j + J

!

!
i,j

Si,j · Si,j+1• Open boundaries in y-direction (multi-leg ladder)• Consider odd L and

J ! ! 1, α limit• Edge spins form triplet | ↑", |0", | ↓" on each rungif

T ! 0.754J ! exp(−L/6.01)• Take first two terms (leg interaction) into accountas weak perturbationEffective Hamiltonian• 9x9 Effective Hamiltonian between two rungs!↑↑ |H | ↑↑# !↑↑ |H | ↑ 0# !↑↑ |H | ↑↓# · · · !↑ 0|H! | ↑↑#Heff = ···!↓↓ |H! | ↑↑# !↓↓ |H! | ↑ 0# !↓↓ |H! | ↑↓# · · ·!

12 non-zero elements among 81

!

!

H =K
!

!
j

!



!↑↑ |H | ↓↓#!↑ 0|H! | ↓↓# ···!↓↓ |H! | ↓↓#S1,j · S2,j

K = 1 or α

• comparing with Hamiltonian of spin-1 dimer J˜ S1 · S2

J˜ = !↑↑ |H! | ↑↑#=K!
j!↑↑ |S1,j · S2,j | ↑↑# = K

!
j

z
z
!↑1 |S1,j| ↑1 #!↑2 |S2,j| ↑2 #z!↑1 |S1,j| ↑1 # : z-component of j-th spin in total Sz=1 sector

QMC Results!
C=!↑ |Sjz | ↑#2
j2

L = 15,17,...,91

1.95
effective interaction

2

1.5

1

1.9
1.85
1.8
1.75

T = 0.02

1.7

0.5

L = 7,9,...,33

00.0001

1.65
1.6

0.001

0.01

0.1

1

0

0.01 0.02 0.03 0.04 0.05 0.06 0.07

T

1/L

No temperature dependence
at T < 0.1

converges to a finite value
for large L

Effective S=1 Chain• Effective S=1 bond alternating chain for J ! ! 1, αHeff =!

J˜i Si · Si+1

J˜ = C or Cα• Quantum phase transition at αc ≈ 0.587• Only spins near surfaces participate

10

Quantum Surface Transition• Phase transitions at finite J’

H-II

1

J’

• No quantum phase transition

in semi-infinite system(i.e. system with only one surface)

AF

0.1
0.08

AF

0.06

H-II

0.04
0.02

H-I

0
0

0.2

0.4

0.6

A

0.8

1

Summary• (Very brief) review on loop cluster quantum

Monte Carlo• Gap estimation by using loop cluster QMC• extended moment method• Topological order in 1D magnets• twist order parameter and its improved

estimator• Quantum surface transition in 2D• perturbation by using QMC

• Randomness effects to spin gapped state

Poster presentations in symposium
by C. Yasuda and M. Matsumoto• Finite-temperature phase transition in quasi-1D

magnets

quantum critical point

Oral presentation in symposium by ST