QMC techniques

Numerical and Analytical Methods for Strongly Correlated Systems
Benasque, Spain, August 24 - September 13, 2014
Anders W Sandvik, Boston University
Quantum Monte Carlo Techniques
(with a focus on quantum spins)
Lecture 1
Stochastic series expansion and ground state projection
Lecture 2
Non-magnetic and critical states in 2D spin systems
Review article on quantum spin systems and
numerical methods: ArXiv:1101.3281
1
Role of numerics/simulations
in studies of many-body ground states and criticality
Obtain definite results for “prototypical” model hamiltonians
(“Ising models of quantum many-body physics”)
- some realized in solid-state materials
- some realizable in cold atoms
- some corresponding to key quantum field theories
- unbiased tests of various analytical calculations
- tools for exploration/discoveries
“Unbiased” methods
(no approximations except finite size)
- exact diagonalization (small systems - be careful!)
- DMRG, 1D systems, recent progress in 2D
- tensor networks, progress in 2D, 3D may be possible
(still convergence issues, can become unbiased in principle)
- QMC, for sign-problem free models, any D, large systems
2
Numerical and Analytical Methods for Strongly Correlated Systems
Benasque, Spain, August 24 - September 13, 2014
Stochastic Series Expansion and
Ground State Projection
Outline
• Path integrals in quantum statistical mechanics
• The series-expansion representation
• Stochastic Series Expansion (SSE) algorithm for the Heisenberg model
• The valence-bond basis for S=1/2 systems
• Ground-state projector algorithm with valence bonds
Reference: AIP Conf. Proc. 1297, 135 (2010); arXiv:1101.3281
Detailed lecture notes on quantum spin models and methods
3
Path integrals in quantum statistical mechanics
We want to compute a thermal expectation value
1
⇤A⌅ = Tr{Ae
Z
H
}
where β=1/T (and possibly T→0). How to deal with the exponential operator?
“Time slicing” of the partition function
Z = Tr{e
H
L
⇤
} = Tr
e
H
l=1
⇥
= /L
Choose a basis and insert complete sets of states;
Z=
0
1
···
L
1
⇥
0 |e
H
|
L
1⇤ · · · ⇥
2 |e
H
|
1 ⇤⇥
1 |e
H
|
0⇤
Use approximation for imaginary time evolution operator. Simplest way
Z⇤
{ }
⌅
Leads to error
0 |1
⇥ H| L 1 ⇧ · · · ⌅ 2 |1
. Limit
⇥ H| 1 ⇧⌅ 1 |1
⇥ H| 0 ⇧
0 can be taken
4
Example: hard-core bosons
H=K=
(a†j ai
Kij =
i,j⇥
+
a†i aj )
ni =
†
ai ai
i,j⇥
Equivalent to S=1/2 XY model
H=
(Six Sjx
2
+
y y
Si Sj )
=
⇥i,j⇤
“World line” representation of
Z⇤
{ }
⌅
0 |1
(Si+ Sj
+
Si Sj+ ),
⇥i,j⇤
⇥ H| L 1 ⇧ · · · ⌅ 2 |1
⇥ H| 1 ⇧⌅ 1 |1
{0, 1}
1
S = ± ⇤ ni = 0, 1
2
z
⇥ H| 0 ⇧
world line moves for
Monte Carlo sampling
Z=
{ }
W ({ }),
W ({ }) =
nK
⇥
nK = number of “jumps”
5
Expectation values
1
⇥A⇤ =
Z
{ }
⇥
0 |e
|
L
1⇤ · · · ⇥
H
2 |e
|
1 ⇤⇥
H
1 |e
A|
0⇤
We want to write this in a form suitable for MC importance sampling
{ }
⇧A⌃ =
A({ })W ({ })
{ }
⇥ ⇧A⌃ = ⇧A({ })⌃W
W ({ })
W ({ })
A({ })
For any quantity diagonal in the
occupation numbers (spin z):
1
A({ }) = A( n ) or A({ }) =
L
=
=
weight
estimator
L 1
A( l )
l=0
Kinetic energy (here full energy). Use
Ke
K
⇧ 1 |Kij | 0 ⌃
1
Kij ({ }) =
⇥ {0,
}
⇧ 1 |1
K| 0 ⌃
K
1
0
1
Average over all slices → count number of kinetic jumps
⇥Kij ⇤ =
⇥nij ⇤
,
⇥K⇤ =
⇥nK ⇤
⇤K⌅ ⇥ N
⇤nK ⌅ ⇥ N
There should be of the order βN “jumps” (regardless of approximation used)
6
Including interactions
For any diagonal interaction V (Trotter, or split-operator, approximation)
e
H
=e
K
e
V
+ O(
Product over all times slices →
W ({ }) =
nK
exp
L
⇤1
l=0
Vl
2
⇥
)⇥⌅
l+1 |e
H
| l⇧
e
Pacc = min
Vl
⇤
⌅
2
l+1 |e
exp
K
Vnew
Vold
| l⇧
⇥
,1
⌅
The continuous time limit
Limit Δτ→0: number of kinetic jumps remains finite, store events only
Special methods (loop
local updates (problem when Δτ→0?)
and worm updates)
•consider probability of inserting/removing
developed for efficient
events within a time window
sampling of the paths ⇐ Evertz, Lana, Marcu (1993), Prokofev et al (1996)
Beard & Wiese (1996)
in the continuum
7
Series expansion representation
Start from the Taylor expansion e
⇥
( ⇥)
Z=
n!
n=0
H
=
⇥
(
n=0
n
{ }n
⇤
)n n
H
n!
(approximation-free
method from the outset)
0 |H| n 1 ⌅ · · · ⇤ 2 |H| 1 ⌅⇤ 1 |H| 0 ⌅
Similar to the path integral; 1
H ⇥ H and weight factor outside
For hard-core bosons the (allowed) path weight is W ({ }n ) = ⇥ n /n!
For any model, the energy is
E
=
=
C = ⇥n2 ⇤
1
Z
this is the operator we “measure”
⇥
( ⇥)n
n!
n=0
1
Z
⇥
{ }n+1
( ⇥)n n
n!
⇥
n=1
⇥n⇤2
⇥n⇤
⇤
{ }n
0 |H| n ⌅ · · · ⇤ 2 |H| 1 ⌅⇤ 1 |H| 0 ⌅
one more “slice” to sum over here
⇤
0 |H| n
⇤n⌅
1 ⌅ · · · ⇤ 2 |H| 1 ⌅⇤ 1 |H| 0 ⌅ =
⇥
relabel terms to “get rid of” extra slice
From this follows: narrow n-distribution with ⇥n⇤
N , ⇥n
N
8
Fixed-length scheme
• n fluctuating → varying size of the configurations
• the expansion can be truncated at some nmax=M (exponentially small error)
• cutt-off at n=M, fill in operator string with unit operators H0=I
=
- conisider all possible locations in the sequence
- overcounting of actual (original) strings, correct by combinatorial factor:
✓
M
n
◆
1
n!(M n)!
=
M!
Here n is the number of Hi, i>0 instances in the sequence of M operators
Z=
X X (
{↵}M {Hi }
)n (M
M!
n)!
h↵0 |Hi(M ) |↵M
1 i · · · h↵1 |Hi(1) |↵0 i
9
Stochastic Series expansion (SSE): S=1/2 Heisenberg model
Write H as a bond sum for arbitrary lattice
2D square lattice
bond and site labels
Nb
H=J
b=1
Si(b) · Sj(b) ,
Diagonal (1) and off-diagonal (2) bond operators
H1,b
=
H2,b
=
z
z
1
S
S
i(b) j(b) ,
4
+
+
1
2 (Si(b) Sj(b) + Si(b) Sj(b) ).
Nb
H=
JNb
H2,b ) +
4
(H1,b
J
b=1
Four non-zero matrix elements
⇤
i(b) ⇥j(b) |H1,b |
⇤⇥i(b)
j(b)
|H1,b | ⇥i(b)
Partition function
Z=
i(b) ⇥j(b) ⌅ =
⇥
⌅⌅
j(b) ⌅
=
⌅
⇥
( 1)n2
n!
n=0
n
Sn
⇥
1
2
1
2
⇤⇥i(b)
⇤
n
⇧1
p=0
j(b)
i(b) ⇥j(b)
Ha(p),b(p)
|H2,b |
i(b) ⇥j(b) ⌅ =
|H2,b | ⇥i(b)
⇤
Index sequence: Sn = [a(0), b(0)], [a(1), b(1)], . . . , [a(n
j(b) ⌅
=
1
2
1
2
n2 = number of a(i)=2
(off-diagonal operators)
in the sequence
1), b(n
1)]
10
For fixed-length scheme (string length = L now)
⇥ L 1
n
⌅⌅
⇧
⇥
(L n)!
n2
Z=
( 1)
Ha(p),b(p)
L!
p=0
S
L
p 1
Propagated states: | (p)⇥
i=0
Ha(i),b(i) | ⇥
⇤
W ( , SL ) =
⇥
2
⇥n
(L
n)!
L!
W>0 (n2 even) for bipartite lattice
Frustration leads to sign problem
In a program:
s(p) = operator-index string
• s(p) = 2*b(p) + a(p)-1
• diagonal; s(p) = even
• off-diagonal; s(p) = off
σ(i) = spin state, i=1,...,N
• only one has to be stored
SSE effectively provides a discrete representation of the time continuum
• computational advantage; only integer operations in sampling
11
Linked vertex storage
The “legs” of a vertex represents
the spin states before (below) and
after (above) an operator has acted
2
3
2
3
2
3
2
3
0
1
0
1
0
1
0
1
X( ) = vertex list
• operator at p→X(v)
v=4p+l, l=0,1,2,3
• links to next and
previous leg
Spin states between operations are redundant; represented by links
• network of linked vertices will be used for loop updates of vertices/operators
12
Monte Carlo sampling scheme
Change the configuration; ( , SL )
Paccept
( , SL )
W ( , SL ) Pselect ( , SL
= min
W ( , SL ) Pselect ( , SL
Diagonal update:
[0, 0]p
[1, b]p
W ( , SL ) =
⇥
, SL )
,1
, SL )
⇥
2
⇥n
(L
n)!
L!
Attempt at p=0,...,L-1. Need to know |α(p)>
• generate by flipping spins when off-diagonal operator
Pselect (a = 0
a = 1) = 1/Nb ,
Pselect (a = 1
a = 0) = 1
W (a = 1)
/2
=
W (a = 0)
L n
(b ⇥ {1, . . . , Nb })
W (a = 0)
L
=
W (a = 1)
n+1
/2
Acceptance probabilities
Paccept ([0, 0] ⇥ [1, b])
=
min
Paccept ([1, b] ⇥ [0, 0])
=
min
Nb
2(L
2(L
n is the current power
• n → n+1 (a=0 → a=1)
• n → n-1 (a=1 → a=0)
⇥
,1
n)
⇥
n + 1)
,1
Nb
13
Off-diagonal updates
Local update
Change the type
of two operators
• constraints
• inefficient
• cannot change
winding
numbers
Operator-loop
update
• Many spins
and operators
can be
changed
simultaneously
• can change
winding
numbers
14
Determination of the cut-off L
• adjust during equilibration
• start with arbitrary (small) n
Keep track of number of operators n
• increase L if n is close to current L
• e.g., L=n+n/3
Example
• 16×16 system, β=16
• evolution of L
• n distribution after
equilibration
• truncation is no
approximation
15
Does it work?
Compare with exact results
• 4×4 exact diagonalization
• Bethe Ansatz; long chains
Susceptibility of the 4×4 lattice
• SSE results from 1010 sweeps
• improved estimator gives smaller
error bars at high T (where the
number of loops is larger)
⇐ Energy for long 1D chains
• SSE results for 106 sweeps
• Bethe Ansatz ground state E/N
• SSE can achieve the ground
state limit (T→0)
16
Numerical and Analytical Methods for Strongly Correlated Systems
Benasque, Spain, August 24 - September 13, 2014
Valence bonds and
Ground State Projection
Anders W Sandvik, Boston University
Outline
• The valence-bond basis for S=1/2 systems
• Ground-state projector algorithm with valence bonds
17
The valence bond basis for S=1/2 spins
Valence-bonds between sublattice A, B sites
Basis states; singlet products
(i, j) = (| ⇥i ⇤j ⌅
⌃
| ⇤i ⇥j ⌅)/ 2
A
N/2
|Vr =
(irb , jrb ),
B
r = 1, . . . (N/2)!
b=1
The valence bond basis is overcomplete and non-orthogonal
• expansion of arbitrary singlet state is not unique
| =
fr |Vr
(all fr positive for non-frustrated system)
r
All valence bond states overlap with each other
Vl |Vr ⇥ = 2N
N/2
N = number of loops in overlap graph
Spin correlations from loop structure
3
xi xj +yi yj
⇤i · S
⇤j |Vr ⌅
⇤Vl |S
(
1)
4
=
0 (i,j in different loops)
⇤Vl |Vr ⌅
(i,j in same loop)
More complicated matrix elements
(e.g., dimer correlations) are also
related to the loop structure
K.S.D. Beach and A.W.S.,
Nucl. Phys. B 750, 142 (2006)
|Vl
|Vr
Vl |Vr ⇥
18
Projector Monte Carlo in the valence-bond basis
Liang, 1991; AWS, Phys. Rev. Lett 95, 207203 (2005)
(-H)n projects out the ground state from an arbitrary state
( H)n | ⇤ = ( H)n
S=1/2 Heisenberg model
H=
i,j⇥
⌅i · S
⌅j =
S
i
ci |i⇤ ⇥ c0 ( E0 )n |0⇤
Hij ,
Hij = ( 14
i,j⇥
⌅i · S
⌅j )
S
Project with string of bond operators
n
⇥
(r = irrelevant)
Hi(p)j(p) | ⇥
r|0⇥
(a,d)
{Hij } p=1
Action of bond operators
Hab |...(a, b)...(c, d)... = |...(a, b)...(c, d)...
1
Hbc |...(a, b)...(c, d)... = |...(c, b)...(a, d)...
2
(a,b)
A
(c,b)
B
(i, j) = (| ⇥i ⇤j ⌅
(c,d)
A
⌃
B
| ⇤i ⇥j ⌅)/ 2
Simple reconfiguration of bonds (or no change; diagonal)
• no minus signs for A→B bond ‘direction’ convetion
• sign problem does appear for frustrated systems
19
Expectation values: A⇥ = 0|A|0⇥
Strings of singlet projectors
n
Pk =
Hik (p)jk (p) ,
k = 1, . . . , Nbn
(Nb = number of interaction bonds)
p=1
We have to project bra and ket states
k
g
Pk |Vr ⇤ =
⇤Vl |Pg =
k
g
Wkr |Vr (k)⇤ ⇥ ( E0 )n c0 |0⇤
⇤Vl (g)|Wgl ⇥ ⇤0|c0 ( E0 )n
6-spin chain example:
A⇥ =
Vl |Pg APk |Vr ⇥
g,k
g,k
=
g,k
Wgl Wkr Vl (g)|A|Vr (k)⇥
g,k
Vl |
A
|Vr
Vl |Pg Pk |Vr ⇥
Wgl Wkr Vl (g)|Vr (k)⇥
- Monte Carlo sampling
of operator strings
- Estimators based on
transition graphs
20
More efficient ground state QMC algorithm → larger lattices
Loop updates in the valence-bond basis
AWS and H. G. Evertz, PRB 2010
Put the spins back in a way compatible with the valence bonds
√
(ai , bi ) = (↑i ↓j − ↓i ↑j )/ 2
and sample in a combined space of spins and bonds
A
|
|
Loop updates similar to those in finite-T methods
(world-line and stochastic series expansion methods)
• good valence-bond trial wave functions can be used
• larger systems accessible
• sample spins, but measure using the valence bonds
21
i =1 2 3 4 5 6 7 8
4
0
9
13
6
0
0
4
13
0
9
14
s(p)
i =1 2
Finite-temperature QMC
(world lines, SSE,...)
tr{e
H
}=
X
n
n
n!
h↵|( H) |↵i
boundaries
by
ANDERS W. open
SANDVIK
AND HANScapped
GERD EVERTZ
valence bonds (2-spin singlets)
[AWS, HG Evertz, 2010]
periodic time boundary conditions
• Computer implementations similar
FIGURE 61. A linked-vertex SSE configuration with one loop tra
“orientations”, along with the corresponding operator-index sequences
flipped, and operators are changed, diagonal ↔ off-diagonal, each ti
change of an operator visited twice). Every vertex leg (spin) belongs
acted upon by any operator (here the one at i = 1) can also be regarde
11
10
9
8
7
6
5
4
3
2
1
0
p
T>0 and T=0 algorithms
of VBs. Expectation values are evaluated at the midpoint indic
ces" are also indicated. There are three loops, part of which co
side-by-side
patible spin states %Zlj&, &Zrj '. The spins at the four operators !v
columns of filled and open circles represent ↑ and ↓ spins of c
Ground state projection
arcs to the left and right indicate VB states %Vl&, &Vr' and the
system with m = 2.
tributing to %%&!−H"2m&%' for a four-site m
FIG. 4. !Color online"↵A VB-spin-operator configuration
n
X
↵
Trial state can conserve relevant
ground state quantum numbers
(S=0, k=0,...)
accomplishes all these things automatically. This class o
duced as a generalization of a cluster algorithm for the Isin
flipped clusters take the form of loops; the classical six-ve
tive world line system for the S = 1/2 Heisenberg model c
Suzuki-Trotter decomposition is exactly equivalent to an
and the loop update for it was therefore at the same time
sical cluster update to a quantum mechanical system. Th
applied also to continuous-time world lines [179] as well a
in the SSE method [190]. The improvements in performa
are enormous (as in the classical case, leading to a much
f f h |( H) |↵i
22
Numerical and Analytical Methods for Strongly Correlated Systems
Benasque, Spain, August 24 - September 13, 2014
Nonmagnetic and Critical Ground States
of 2D Quantum Spin Systems
SA
SB
ective description of the rotationally invariant Néel vector mp
s in terms of two large
#iThere
"j )/ is2an effective antiferrosponding to the sum of the spins
sublattices.
~ion the two ("
i #j
S
between these spins, leading to a singlet ground state and a “tower” of quantum rotor
spin S = 1, 2, . . . at energies ΔS ∼ S(S + 1)/N above the ground state.
23
Outline
Conventional quantum phase transition in 2D antiferromagnets
- Néel to non-degenerate quantum paramagnet
Unconventional transition (deconfined quantum criticality?)
- Néel to valence-bond-solid (4-fold degenerate ground state)
- Sign-free QMC realization: “J-Q” models
SA
SB
Effective description of the rotationally invariant Néel vector ms in terms of two large
Studies
J-Q
models
rresponding
to the of
sumcriticality
of the spins on in
the two
sublattices.
There is an effective antiferrong between these spins, leading to a singlet ground state and a “tower” of quantum rotor
- Finite-size scaling
al spin S = 1, 2, . . . at energies ΔS ∼ S(S + 1)/N above the ground state.
Universality: Correspondence in frustrated spin models?
- comparisons
with
recent
results
forspin-wave
J1-J2 theory.
Heisenberg
model
for the
symmetry breaking
that is
the starting
point for
In
namic limit, the direction of the ordering vector is fixed (as the time scale
th its rotations diverges [169]), and the quantum rotor-states are then in
ccessed. They are neglected in standard spin-wave calculations (discussed
om the outset because the order is by construction locked to the z direction.
24
Starting point: S=1/2 antiferromagnetic Heisenberg model
H = J
i,j⇥
Si · Sj
Sublattice magnetization
1
m
⌃s=
N
N
⌃
i Si ,
i
= ( 1)xi +yi (2D square lattice)
i=1
Long-range order: <ms2> > 0 for N→∞
Quantum Monte Carlo
- finite-size calculations
- no approximations
- extrapolation to infinite size
L⨉L lattices up to 256⨉256, T=0
Reger & Young 1988
ms = 0.30(2)
60 % of classical value
AWS & HG Evertz 2010
ms = 0.30743(1)
25
T=0 Néel-paramagnetic quantum phase transition
Example: Dimerized S=1/2 Heisenberg models
• every spin belongs to a dimer (strongly-coupled pair)
• many possibilities, e.g., bilayer, dimerized single layer
strong interactions
weak interactions
Singlet formation on strong bonds ➙ Néel - disordered transition
Ground state (T=0) phases
s
= spin gap
3D classical Heisenberg (O3) universality class; QMC confirmed
Experimental realization (3D coupled-dimer system): TlCuCl3
26
L=8
L = 16
L = 32
(a)
g = 1.0
g = 1.5
g = 1.9
g = 3.0
0.25
(b)
(a)
L=8
L = 16
L = 32
0.10
0.05
1.5
2.0
g
0.20
g = 1.0
g = 1.5
g = 1.9
0.10g = 3.0
0.15
0.05
0.10
3.0
0.00
0.25
0.05
0.00
1.0
0.20
2.5
<ms2>
0.10
<ms2>
<ms2>
0.15
<ms2>
0.15 of finite-size scaling scaling studies; dimerized Heisenberg
Example
0.05
0.15
0
(b)
0.05
0.1
0.15
0.2
0.25
1/L
ρsL
FIGURE 5. QMC results for the squared sublattice magnetization
in theto
two-dimensional
Heisenberg
g
=
J
/J
According
theory, spin
stiffness
2
1
0.00
0.00
L = results
8
model1.0
with 1.2
columnar
dimerization.
g for0.25
different lattice sizes
2.0
3.0
0 versus
0.1
1.5
2.5 (a) shows
0.05the coupling
0.15ratio0.2
g
the1/L
critical
should
scale
and (b) shows the size dependence
for several
A quantum
phase point
transition
where the
Néel
L = 16values of g. at
order vanishes occurs at g ≈ 1.9.
according
to (T=0)
L = magnetization
32
GURE 5. QMC1.1
results for the squared sublattice
in the two-dimensional
Heisenberg
L=
64 the coupling ratio g for different lattice sizes
del with columnar dimerization. (a) shows results
versus
1
According
to
theory, spin stiffness
128g. A quantum ⇢
(b) shows
phase⇠
transition
where
theconstant
Néel
1.2the size dependence for several
L =values
8 L = of
!
L⇢
s
s
1.0
thetheory—the
criticalLpoint
should
scale in
er vanishes
occurs at g ≈ 1.9. treatments
σ -model
renormalization-group
of one suchatfield
nonlinear
L = 16
to (T=0)
2+1 1.1
dimensions [5, 84]. BasedL =on32symmetry according
arguments alone,
one would then expect
Allows
accurate
determination
0.9 to be in the universality class of the 3D
the transition
classical
Heisenberg
model. Thereof the
L = 64
critical
point
crossings)
σ -model
in
ormalization-group
treatments
ofthe
one
such field theory—the
nonlinear
are, however, subtle
issues in
mapping,
and(curve
QMC
simulations
are
L =quantum-classical
128
s
s
1.0
1 dimensions
[5, 84]. Based on symmetry arguments alone, one would then expect
0.8
ρsL
1
⇢ ⇠ ! L⇢ constant
We will see L
examples of such comparisons
1.88 to test
1.90 various
1.92predictions.
1.94
1.95therefore needed
transition
be in the
universality
classfield
of the
3D classical
model.
g and
betweentoresults
of simulations
theories
in Sec.Heisenberg
5. While results
forThere
the transition
Allows
accurate
of the
, however,
subtle issues
in the
quantum-classical
QMC simulations
are several
0.9 bilayer
in the
(a) [85]
and
columnar dimermapping,
(b) [86]and
systems
indetermination
Fig. 4 (and
drefore
spin stiffness
(in
the
x various
direction)predictions.
multiplied byWe
thewill
system
needed
to
test
see
examples
of (curve
such comparisons
critical
point
crossings)
other
cases
[87,
88])
are
in
good
agreement
with
the
expectations,
recent studies of
model. The crossing points of these curves for different L tend
ween
results
simulations
field
theories indeviations
Sec. 5. While
for still
the transition
staggered
dimers
(c)and
show
[89]results
that are
not understood.
0.8barsof
io g.the
Error
are much
smaller
than
theunexpected
symbols.
1.88(a) [85]
1.90and columnar
1.92
1.94 (b) [86] systems in Fig. 4 (and several
the bilayer
dimer
g
her cases [87, 88]) are in good agreement with the expectations, recent studies of
staggered
dimers
(c) show
unexpected
deviations
tiffness
(in the
x direction)
multiplied
by the
system [89] that are still not understood.
27
More complex non-magnetic states; systems with 1 spin per unit cell
H = J
i,j⇥
Si · Sj + g ⇥ · · ·
• non-trivial non-magnetic ground states are possible, e.g.,
➡ resonating valence-bond (RVB) spin liquid
➡ valence-bond solid (VBS)
Non-magnetic states often have natural descriptions with valence bonds
i
j
= (⇥i ⇤j
⌅
⇤i ⇥j )/ 2
The basis including bonds of all lengths
is overcomplete in the singlet sector
Spontaneous symmetry breaking
(different from dimerized Hamiltonian)
• non-magnetic states dominated
by short bonds
28
Non-magnetic states from frustrated spin interactions
Quantum phase transitions as some coupling (ratio) is varied
• J1-J2 Heisenberg model is the prototypical example
H=
i,j⇥
= J1
⌅i · S
⌅j
Jij S
= J2
g = J2 /J1
• Ground states for small and large g are well understood
‣ Standard Néel order up to g≈0.45; collinear magnetic order for g>0.6
0
g < 0.45
0.45
g < 0.6
g > 0.6
• A non-magnetic state exists between the magnetic phases
‣ Most likely a VBS (what kind? Columnar or plaquette?)
‣ Some recent calculations suggest spin liquid (but I doubt it...).
• 2D frustrated models are challenging: QMC sign problems
29
VBS states from multi-spin interactions
(AWS, PRL 2007)
The Heisenberg interaction is equivalent to a singlet-projector
Cij =
Cij |
1
4
s
ij ⇥
⇤i · S
⇤j
S
=|
s
ij ⇥,
Cij |
tm
ij ⇥
= 0 (m =
1, 0, 1)
• we can construct models with products of singlet projectors
• no frustration in the conventional sense (QMC can be used)
• correlated singlet projection reduces the antiferromagnetic order
+ all translations
and rotations
The “J-Q” model with two projectors is
H=
J
Cij
ij⇥
Q
Cij Ckl
ijkl⇥
• Has Néel-VBS transition, appears to be continuous
• Not a realistic microscopic model for materials
• Intended to study VBS and Néel-VBS transition (universal physics)
30
fieldindex).
that conveniently
describes
co
6 theglo
from
these
(where
µ
=
x,
⌧
is
a
three-component
space-time
The
field
be expressed
universally as functions ofexcitations
the such
spin
systems
,
and
a
carrying
S
=
1
of
the antiferromagnet.
As w
order
parameter.
This
means
that
the
two
order-p
VBS
states
and
“deconfined”
quantum
criticality
gauge
invarianc
A
is
unrelated
to
the
electromagnetic
field,
but
is
an
internal
µ
kquite
ground
state
of
H
B T /⇢ s .
0 also ha
N´
e
el
state,
expressing
the
spin-wave
fl
distinct
fields
will
have
long-ranged
“statistical”
interactio
the
physical
ob
field
that
conveniently
describes
the
couplings
between
the
spin
Sachdev
Vishwanath,
Balents, Sachdev,
Fisher
(2004)
sis, Read,
it is useful
to(1989),....,Senthil,
have an alternative
description
pattern
of
spin
polarizatio
and
A
is
a
matter
of
choice,
and
the
a
µ
transformation
excitations
the antiferromagnet.
As
we have noted there
above,will
in thebe no
connectedof to
each
other.
Consequently
local
theo
field
8
can
serve
us
equally
well.
The
d
states
above
the
N´
e
el
ordered
state.
For
the
equation
(1)
by
H
=
J
S
·
S
+
g
⇥
·
·
·
i
jthe spin-wave fluctuations in terms of z ↵
N´of
eel more
state, expressing
s
exincludes
only
thea two
order-parameter
(butou
approaches
appears whenfields
we move
ernative
description
is,
in
a
sense,
purely
i,j⇥
and A µ is a matter
of choice, and the above theory for the vector
nise:
principle:
in
quantum critical
points
into
other
pha
fields).
It
is
these
difficulties
that
force
the
necessit
it
does
not
alter
any
of
the
low-energy
h
S
N
fieldNeel-VBS
8 can serve
us equally
well. The distinction between the two
j!i =
1c
!
transition
in
2D
!Si⌘(r,
· Sj " ⌧)A
inwhich
some ofisthese
phases,where
the= emergent
ons
resonate
alternate
description
conveniently
provide
approaches
appears
when
we move out of theory
the
N´eel state
across
and
an
identical
low-temperature
•yields
generically
continuous
optional,
but an
essentialimportant
characterizatio
role
resulting
state
quantum
critical
points
into
other
phases
(as
we
will
see
later):
spinon
ofkey
freedom.
hen•expressed
in terms
of rule”
kBdegrees
T /⇢s . The
step
r i isforthe
position
of
violating the
“Landau
the
phase.where
As we did
S
,
we
can
writ
8
theory of the
in
some
of
these
phases,
the
emergent
A
gauge
field
is
no
longer
µ complex
ractional
spin
ctor
field
8
in termstransition
of
an spinon
S = 1/2fields
The
defined
in constraints
Eq.equation
(1.5)
have
forzz"↵ and
A µ by the
symm
ordering
pattern
inofFig.
1b
stating
1st-order
(6)
ca
characterization of thewhich
‘quantum
order’ of
now
yields
eoptional,
↵ Description
="# but
by anofessential
most
important
diVerence
“gauge”
redundancy.
Specifically
the
local
ro
or
w
. We will
critical
point
with
spinor
field
↵phase
the phase. As we did for S8 , we can write the quantum field Ztheory 
olution
ofA µthe
two orthogonal
vectors
Section
IID. N2
complex
vector)
2
for (2-component
z ↵ and
by
the
constraints
of
symmetry
and
gauge
invariance,
i$$r,d%% r d⌧ |(@ µ iA µ )z ↵ |
⇤
S
z =
gauge
redundancy:
z of ground states
z
→
e
8 =yields
z↵ ↵ z
(3) manifold
urwhich
primary
exnow
B. COUPLED-DIMER
A

Z
an eVective
action
for
N
1
,
nd a valence
1
leaves
the
Néel
invariant
and
hence
is
gaug
2a mode
2
2 vector
2
2 2
This
spin
(✏
@
A
)
.
+
2ions
Paulibetween
that
this
Smatrices.
d rNote
d⌧ |(@
iA µ )zmapping
+ u(|z ↵ |the
) hamiltonian.
µ⌫l ⌫ l Minimiz
z =
µ
↵ | + s|z ↵ |from
2
2e0 time
%
is
the
imaginary
of
freedom.
Here
has coordina
N21 = N22 fi
.eeWeSec.
can VIII).
make a space-time-dependent change ordered stateantiferromagne
dashed
lines
to
1 brevity,
For
we
have
now
used
a
‘relativi
the
spinons
are
coupled
to
a
U$1%
gauge
field
a
1
•
CP
action
(non-compact)
Afield
is a U(1)
symmetric
gauge
&$
in
the
hamiltonian,
but
the
✓(x,
⌧)
2 field
wn to be ob(5)spin-wave
+ 2 (✏ µ⌫l @ ⌫ Al ) .
square-lattice
m;
and
scaled
away
the
velocity
v
& ,analogue
' , . . . tocan
represent
t
of
the
spino
2e0 will use the Greek indicesthe
be underst
- proposed as
i✓ critical theory separating Neel and VBS states
). Our
central
thesis—subs
space-time
indices xnature
,invariant
y(4)
, %physics
introduce
another
spinor
z↵ !
ehave
z ↵ now
VOL
4 MARCH
2008
www.nature.com
For
brevity,
we
used
a
‘relativistically’
notation,
exchange
intera
N-1
ce bond
solid
• SU(N) generalization:
large-N
calculations
for CP
theory
by
a
variety
of
arguments
to
follow—is
crit
and scaled away the spin-wave velocity v ; the values of the couplings ofthat
J /g the
. A num
aks
spin
rota[can be
carried out with similar QMC as SU(2) models]
nged. All physical theory
properties
therefore transition is just theNsim
for must
the Néel-VBS
1+
attice
translaVOL 4 MARCH 2008 www.nature.com/naturephysics
nature physics
31
T=0 Néel-VBS transition in the J-Q model
Ground-state projector QMC calculations
(Sandvik, 2007; Lou, Sandvik, Kawashima, 2009)
Néel order parameter (staggered magnetization)
1
⌅
M=
N
⌅i
( 1)xi +yi S
i
VBS vector order parameter (Dx,Dy) (x and y lattice orientations)
1
Dx =
N
N
( 1) Si · Si+ˆx ,
xi
i=1
1
Dy =
N
N
i=1
( 1)yi Si · Si+ˆy
No symmetry-breaking in simulations; study the squares
⌅ ·M
⌅ ⇤,
M 2 = ⇥M
D2 = ⇥Dx2 + Dy2 ⇤
Finite-size scaling: a critical squared order parameter (A) scales as
A(L, q) = L
(1+ )
f [(q
qc )L1/⇥ ]
Data “collapse” for different system
sizes L of AL1+η graphed vs (q-qc)L1/ν
coupling ratio
Q
q=
J +Q
32
J-Q2 model; qc=0.961(1)
s
d
= 0.35(2)
= 0.20(2)
J
Q2
⇥ = 0.67(1)
J-Q3 model; qc=0.600(3)
s
d
= 0.33(2)
= 0.20(2)
⇥ = 0.69(2)
Exponents universal
(within error bars)
Comparable results for
honeycomb J-Q model
Alet & Damle, PRB 2013
Block, Melko, Kaul, PRL 2013
J
Q3
Exponents drift for large L
Kawashima et al, PRB 2013
- weak first-order transition?
- or large scaling corrections?
33
Universality of J-Q physics: Frustrated spin models
H=
i,j⇥
= J1
⌅i · S
⌅j
Jij S
= J2
g = J2 /J1
Until recently, most calculations indicated VBS around J2/J1=1/2
0
g < 0.45
0.45
g < 0.6
g > 0.6
Recent DMRG calculations
claim a spin liquid
1) Jiang, Yao, Balents (PRB 2012)
2) Gong, Zhu, Sheng, Motrunich,
Fisher (arXiv 2014)
34
ns and PVB textures whose decay length grows strongly with
ge PVB order in the two-dimensional limit. The dimer-dimer
B order.
For 0.44 <ordered
J2 /J1 < 0.5,
both spin
order,
dimer
Plaquette
phase
and
quantum
s and appear to scale to zero with increasing system width,
near-critical behavior. We compare and contrast
our results
1
1
1
spin liquid in the spin- 12 J1 -J2 square Heisenberg model
Shou-Shu Gong , Wei Zhu , D. N. Sheng1 , Olexei I. Motrunich2 , Matthew P. A. Fisher3
Department of Physics and Astronomy, California State University, Northridge, California 91330, USA
2
Department of Physics, California Institute of Technology, Pasadena, California 91125, USA
arXiv:1311.5962v1
3
Department of Physics, University of California, Santa Barbara, California 93106-9530, USA
2
ms
ms
2
energies. On1the other hand, recent DMRG
studies[45–47] of
0.1
te
0.1
We study
the spin- 2 Heisenberg model on the
square
lattice with first- and
second-neighbor antiferromagnetic
0.2
(a)
another
bipartite
frustrated system—the J1 -J2 spin-1/2 hon0.08


iJ1 and
J
and
might
ms eycomb
2 , which possesses a nonmagnetic region that has been debated for many years
0.4 interactions
0.06 2
lattice Heisenberg model—found a PVB phase in the
m s, with
ro
realize thenonmagnetic
interesting region,
Z2 spin
liquid
(SL).
We
use
the
density
matrix
renormalization
group
approach
0.04
with a possible SL phase between the
0.15
0.075
0.02 cylinders
in
explicit
of SU (2)
spinarotation
symmetry
and study the model accurately on open
⇥ implementation
N´eel and PVB phases[47]
or with
direct N´eel
to PVB transi0
0.3 with different
rconditions.
Withquantum
increasing
J2point[45–
, we find a N´eel phase,0 a 0.1
plaquette
(PVB)
tion boundary
characterized
by a deconfined
critical
0.2 0.3 valence-bond
0.4 0.5

J
These
studies[46,
also foundspin
that liquid
in the nonmagnetic
2
phase with49].
a finite
gap, and47]
a possible
in a small region of J0.1
two phases. From the
Tc
spin
2 between these
⇥
0.05
region
theofconvergence
of DMRG
in
wider systems,
which is that the N´eel order vanishes at J2 /J1 ' 0.44. JJ =0.35
finite-size
scaling
the
magnetic
order
parameter,
we
estimate


x=0.40
0.2
controlled
by
the
number
of
states
kept,
is
crucial
for
deterJ =0.44
For 0.5 < J2 /J1 < 0.61, we find dimer correlations and PVB textures whose
0.05 decay length grows strongly with J =0.46
imining
the
true
nature
of
the
ground
state.
increasing system width,gapless
consistent with a long-range PVB order in the two-dimensional limit. The dimer-dimer J =0.50
fiIn
this
Letter,
we reexamine the J0.025
2 square lattice
0.1 correlations
reveal the s-wave
character of the1 -J
PVB
order. For 0.44 < J02 /J1 < 0.5, both spin order, dimer J =0.55
region
NeelHeisenberg
3–
model for J2 < PVB
0.61 using DMRG algorithm
0zero with
0.1 system
0.2
0.05increasing
0.15 width,
0.25
order,
and
spin
gap
are
small
on
finite-size
systems
and
appear
to
scale
to
iwith explicit implementation of the SU (2) spin rotation
1/L our results
which
is
consistent
with
a
possible
gapless
SL
or
a
near-critical
behavior.
We
compare
and
contrast
symmetry[50]
(we
do
not
study
the
well
known
stripe
antil0 with earlier numerical studies.
0
0.2
0.3
0.4
0.5
ferromagnetic
phase
at
larger
J
).
We
find
accurate results
2
al
(b)
on cylindersJ with system width up to 12 ⇠ 14 lattice spac2
st
PACS numbers:
ings 73.43.Nq,
by keeping75.10.Jm,
as many 75.10.Kt
as 36000 optimal U (1)-equivalent
0.1
rstates. We find N´eel phase
below J2 ' 0.44 and a nonmag1
FIG. 1: (color online) Phase
diagram
spinsquare
up
Order-parameter
extrapolations
indicate
netic region
forof0.44
<2 JJ21 -J
< 20.61
by Heisenfinite-size
scaling of
berg
model
for
J
<
0.61
obtained
by
our
SU
(2)
DMRG
studies.
2
r- Introduction.—Quantum
the magnetic
order parameter.
In the
nonmagnetic
spin liquid
(SL) is an
exotic
state region, we
0.1
disordered
but
extrapolations
are
With growing J2 , the phase,
model has
a N´
eorder
el phase
J < 0.44 and a
establish
a PVB
for
Jfor


2 > 20.5—in contrast to the previonof matter
where a spin system does not form magneti2
2
2
2
2
ms
2
2
PVB phase for 0.5 < ous
J2 proposal
< 0.61. [37]
Between
two observing
phases, there
0.4
of a Zthese
that the PVB de2 SL—by
callyisordered
state or
break
lattice
even
at system
zero width. We
a small gapless
region
exhibits
no orderwith
in our
calculations,
cay that
length
growssymmetries
strongly
increasing
istemperature[1].
consistent with Understanding
a gapless
SL. The
mainorder
panel
N´eelplaquette[30]
order paidentify
the PVB
as shows
theis
s-wave
spin
liquids
important
in by studyrametermagnetic
ms and spin
gap
in the
thermodynamic
limit.
ing dimer-dimer
correlations.
For
0.44
<
Jinset
0.3
Tand
2 < 0.5, we find
frustrated
systems
may
also hold clues
toThe
underis a sketch for an RC4-6
Jpin shows
the modified crystal
odd verthatcylinder;
the magnetic
order, valence-bond
(VBC) orders,
standtical
non-Fermi
liquid
of
doped
Mott
materials
and
high-T
c
as well
as spin excitation
vanish
with increasing
sysNN bonds providing
the boundary
pinninggap
for all
dimer
orders.
1)
ms
ms
no reliable close to a critical point
- most likely dqc point as in J-Q
0.075
⇥

superconductivity of strongly
correlated
systems[2].
The SL
ex-in agreement
tem width,
which suggests
a possible gapless
0.2
the VMC
or a near-critical
behavior.
citing
propertiesresults
of spinwith
liquids
suchresults[41]
as deconfined
quasipartiDMRG
should
be
compared
incaldetail
We
consider
both
torus
and
cylinder
samples
in
DMRG
have been
revealed
in many
less fractional
SL[27, 28, statistics
40, 41].
However,
true nature
of the artifiquant-cles and
culations,
but all the
the phases
are established
based on high ac0.1
cially
constructed
systems
such
as
quantum
dimer
models[3–
tum
phase
has
remained
unresolved.
curacy results on cylinders[51]. We use two cylinder geomes,
modeltries.
in the
limit[6–10],
and
Ki- with closed
Theeasy
first isaxis
thethe
rectangular
cylinder
(RC)
ng5], Kagome
Recentspin
large-scale
DMRG
study
of
J1 -J2 square
lattice
boundary inofthefinding
y direction
open boundaries
taev model[11]. The possibility
spinandliquids
in real- in the x di-

L
10
0.05 ⇥

FIG. 2: (color online) (a) m2s plotted vs 1/L for RCL-2L cylinder
with L = 4, 6, 8, 10, 12, 14; lines are polynomial fits up to fourth
order. The inset is J2 dependence of the obtained magnetic order
2
in the 2D limit m2s,1 . (b) Same data shown as log-log plot of m0.025
s
versus width L.
with J-Q QMC
gapless
region PVB
Neel
35
... New version of the paper [PRL 113, 027201 (2014)]
Plaquette ordered phase and quantum phase diagram in
the S=1/2 J2-J1 square Heisenberg model
“The critical exponents obtained from the finite-size spin and dimer
correlations could be compatible with the deconfined criticality”
Conclusion from studies of J-Q and frustrated square lattice
- the J-Q model can mimic the behavior of (some) frustrated systems!
- many more insights into deconfined criticality and VBS states obtained
by large-scale QMC studies of J-Q models
K. Harada et al, PRB 88, 220408 (2013)
M. Block, R. Melko, R. Kaul, PRL 111, 137202 (2013)
S. Pujari, K. Damle, F. Alet PRL 111, 087203 (2013)
Y. Tang, AWS, PRL 110, 217213 (2013)
S. Jin, AWS, PRB 87, 108040 (2013)
AWS, PRB 84, 134407 (2012)
A. Banerjee, K. Damle, F. Alet, PRB 83, 235111 (2011)
.....
36