Computational Approaches toQuantum Critical Phenomena(2006.7.17-8.11) ISSPFermion Simulations
July 31, 2006
Univ. Tokyo

M. Imada
ISSP, Kashiwa

collaboration
T. Kashima,
Y. Noda,H. Morita,
T. Mizusaki,
S. Watanabe
I.V. Solovyev
Y. ImaiY. Otsuka

M. IMADA

Numerical Methods for Lattice Fermions
exact diagonalization
high-temperature expansion
DMRG
variational Monte Carlocellular (cluster) DMFT
quantum Monte Carlo
auxiliary fieldworld line
PIRG
Gaussian basis

small clusters
high to moderate T
1D systems
trial wave functions
small cluster + DMFnegative sign

PIRG: nonperturbative but systematic approach from U=0

the only available approach for frustrated models,M. IMADA
complex systems with orbital at T=0 for large
sizePath-Integral Representation

Φ = lim[exp[ −∆ β H ]] Φ 0pp→∞

Φ' = exp[−∆ β H ] ΦPractically,taking small ∆β,but ∆βp large≈ exp[−∆ β H 0 / 2]exp[−∆ β H1 ]exp[−∆ β H 0 / 2] Φ + O( ∆ β 3 )Hubbard model

H = H 0+H 1
H 0 = Σσ tij c†iσ cjσ + H.cij

H 1 = U Σｎi↓ｎi↑example’tt M. IMADA
UM. IMADA

Matrix representation of Slater determinant
M Fermions on N-sites lattice

inner product|Φ〉

MN(Φσ )im

Example; plane wave state

M. IMADASDW Hartree-Fock State

M. IMADAOperation of Kinetic Energy Projection
Φ = exp[−∆ β H 0 ] Φ'bilinear tight-binding form

H 0 = Σ [tij c†iσ cjσ + H.c + δijµ c†iσ cjσ ]Matrix representationof exp[−∆ β H 0 ] ⇒ b0NN

In Slater determinant
representation, exponentialof kinetic energy generates
another Slater detereminantsMΦΦ′Φ
M. IMADAOperation of Coulomb Energy Projection
Φ ' = exp[−∆ β H1 ] Φ ??

H1†
= UΣ ci↑ ci ↑†ci↓ ci ↓Projection by the interaction term
is transformed to a sum oftwo Slater determinantsby Stratonovich-Hubbard transformation
M. IMADA

Stratonovich-Hubbard Transformation

Discrete transformation
Auxiliary Ising variable sdecompositionn↑ n↓ ⇒ n↑ × n↓
M. IMADA

Stratonovich-Hubbard transformation
Stratonovich-Hubbard transformation

Φ = exp[−∆ β H1 ] Φ = ∑ Φ( s )'s＋branching occurs
M. IMADA

Basis Generation in Slater Determinant Representation

Φ = lim[exp[−τH]]
Φ0p →∞pexp[−τ H ] ≈ exp[−τ H K ]exp[−τ H U ]In Slater determinant representation
Kinetic energyΦ′Φ = exp[−τ H K ] Φ'ΦInteraction energy
Stratonovich-Hubbard transformation

Φ ' = exp[−τ H U ] Φ = ∑ Φ( s ) ; branchings＋M. IMADA

Ground State Algorithm
by Auxiliary Fields of Path Integral
ground stateΦ = lim[exp[ −∆ β H ]] Φ 0pp→∞ΦΦ0e0−τHM. IMADA

Matrix Elements in Canonical Ensemble
Trace out Fermion degrees of freedomleaving auxiliary fields

density matrix

M. IMADAPath integral in grand canonical ensemble
valid at finite temperature

β = p∆ β

M. IMADA

Full summation over path integral

summation over all {s}M. IMADA

Numerical Algorithms
Summation over all
the Stratonovich-Auxiliary variablerequires 2Npterms summations

Monte Carlo sampling

quantum Monte CarloRenormalization(projection andtruncation) ofgenerated basis

path-integral
renormalizationgroup
M. IMADA

negative sign problem

Anticommuting Fermions

W(A) = 〈i | e −τ H | j〉〈 j | e −τ H | k〉〈k | eii−τ H....A....| i〉W(1) < 0

Average sign decreases exponentially withincreasing Ν, β∑ W (A) 0A =≈∑ W (1) 0ss

Large statistical errorss

Difficulty in QMC
M. IMADAPath-IntegralRenormalization Group(PIRG)
Kashima et al. JPSJ 69 (2000)2723; 70(2001)2287

M. IMADA

Numerical Framework of PIRG
Kashima et al. JPSJ 69 (2000)2723; 70(2001)2287

Optimize Φ at fixed LLΦ = ∑ ci ϕ iIncrease LiExtrapolationE=Φ H ΦΦΦA =ϕi A ϕ j , ϕi ϕ j
ex. Slater determinantΦ AΦΦΦ

variationalEasily computable
No negative-sign problem
M. IMADA

Basis Generation in Slater Determinant Representation

Φ = lim[exp[−τH]]Φ
0p →∞ppath integral

exp[−τ H ] ≈ exp[−τ H t ]exp[−τ H U ]In Slater determinant representation
Kinetic energyΦ′Φ = exp[−τ H t ] Φ'ΦInteraction energy
Stratonovich-Hubbard transformation

Φ ' = exp[−τ H U ] Φ = ∑ Φ( s ) ; branchings＋nonorthogonal basis

M. IMADAWF Renormalizationin the direction to imaginary timeΦ = lim[exp[−τ H ]] p Φ 0p→∞

exp [− τH ] Φ 0p

Projection to the ground state

PROJECTION
Select the subspace which lowers the energy

orTRUNCATION

M. IMADAComparison with Orthogonal Basis
Hartree-Fock energy

with
Gutzwiller factor(VMC)

Orthogonal basis= Configuration Interaction

Optimized non-orthogonal basis
M. IMADARenormalization procedure
At a fixed dimension, L in restricted Hilbert spacein some representation, find the “ fixed point ” :ϕ i ϕ j and ϕ i H ϕ j“effective Hamiltonian”

L∑projection &Φ = ci ϕ i ground state oftruncation“effective Hamiltonian”
iw.f. renormalization,
hamiltonian matrix renormalization
Renormalization in energy; filtering out high energy modec.f. in usual RG, one finds a fixed point rep. in thescaled but essentially fixed rep.
M. IMADAExtrapolation to L→∞∆E =(E2 − E2)

E − E0 ∝ ∆ E/ E

2: Extrapolation with energy variance: L large

6 × 6 lattice at half filling

VMC

S=1n=0S=0E = −17.693
1En =exact = −17.6947

2

M. IMADAQuantum number projection
variational ground states in the restricted Hilbert space
do not necessarily preserve the original symmetries of H
How to restore the symmetries ?

Quantum number projection operator

Mizusaki and Imada
PRB69, 125110 (2004)
M. IMADAMomentum Projection

K; momentum operator
j spatial translation

M. IMADA

Spin rotation
Wigner’s D function
Spin rotation
Euler angle

N0 = Sz
M. IMADA

Spin projection

Legendre polynomial
M. IMADAQuantum-Number Projected PIRG (QP-PIRG)
Mizusaki and Imada
PRB69, 125110 (2004)

is replaced with

Yrast state;lowest energy state with specified quantum number
M. IMADAQuantum Number Projection + PIRG
Spin, Momentum, …
QMC -1.8574(14)QP-PIRG -1.85790（2）PIRG+QP
QP-PIRG

VMC

S=1
S=0

QMC (Assaad)

Excited States
Energy Dispersion

Yrast states
M. IMADAExcitation Spectra, Dispersion
-33.2

-33.2

S=0
S=0 exact

-33.4

S=1
S=1 exact-33.4

EE-33.6
-33.6

-33.8
-33.8

-34
-34.2

(0,0) (π/2,0) (π,0) (π,π/2) (π,π) (π/2,π/2)k-34(0,0) (π/2,0) (π,0) (π,π/2) (π,π) (π/2,π/2)k

4×4 Hubbard
U=5.7, t=1, t’=0.5
M. IMADA

Grand Canonical Ensemble
Watanabe, MI;
JPSJ 73 (2004) 1251

H = H − µ Nˆ eck ↑ck−k ↓+k{cParticle-hole transformation
Transformed Hamiltonian:d⎛U⎞
H = H t + HU − ⎜ + µ ⎟ N
⎝4
⎠

()

⎛UH t = − ∑ tij c c j + c c + ⎜ − µ ⎞⎟∑ ci+ ci⎝2⎠ iij⎛U⎞+ ∑ tij (d i+ d j + d +j d i ) + ⎜ + µ ⎟∑ d i+ d i⎠ i⎝2ij+i+j i

H U = −U ∑ c c d d i+i ii+i

M. IMADA

Extended basis⎛+⎞~φ = ∏ ⎜ ∑ [φ ]ik ci ⎟ 0k =1 ⎝ i =1
⎠M

2Ncanonicalc~i + =ci+{d+i− N

for i = 1,..., N
for i = N + 1,...,2 N

grandcanonicalck+ d k = ck+↑ c−+k ↓[φ ]ik :cf. superconductivity
Yokoyama & Shiba(1988)

Total electron number

N e = ∑ ci+σ ciσ =N + ∑ ci+ ci − d i+ d iiσiM. IMADA

Several Applications

M. IMADAHubbard Model

M. IMADAPhase Diagram of Mott
Transition in the 2D Hubbard model; T=0PIRG numerical results
t = 1, t ' = −0.2

8Mott Insulator

U /t

6filling control
TMO

42N

4× 4
6×6
8×8
10 × 10

Filling-Control

organics

Metalbandwidth
Metalcontrol

1st order transition

V shape boundary
0Reduction-2 -1.5 -1 -0.5 0 0.5
1 1.5 2 of jump for larger t’
Watanabe, MIµµ/ t

M. IMADAJPSJ73(2004)1251

Mott transition in experiments
κ-(ET)2Cu[N(CN)2]Cl

V2O3chargeorderRNiO
3(DI-DCNQI)2Ag
PRL93(2004)216408

IM
Limelette et al.
Kanoda et al.

Reduction of Tc
for larger frustration
pyrochlore; continuous
M. IMADAPhase Diagram of Frustrated
Hubbard Model at ½ Filling

Frustrated
square lattice
Kashima et al.
JPSJ 70 (2001) 3052

Spin melting
before charge melting

Quantum melting
Anisotropic by charge fluctuations
HF;Vollhardt et al.

triangularlatticet’

Morita et al.

t

JPSJ 71(2002) 2109

UIMADA
Nonmagnetic Mott insulatorM.phasespin liquids near MIT

No symmetry
PM breaking⇔ dimer phase(Sachdev, Sushkov,...)

Mizusaki & MI,
PRB74, 014421 (2006)

M. IMADANature of spin liquid phase: Spin-gapless phase
gapless spin excitations∆E ~ 1/N
Triplet excitation gap collapses
with increasing system size

Finite uniform susceptibility

Mizusaki & MI,
PRB74, 014421 (2006)

M. IMADAk-dependence of gapless excitation
width of excitation energy

4E
3

k1/N

2
1
S=0

Gapless excitations withdouble hierarchy structure
cf. kagome, pyrochlore;
Chalker et al., Canals, Lacrois
Harris et al., Tsunetsugu

Very small dispersion (incoherent excitation)?
orcontinuum of excitation (spinon Stoner excitation)?
M. IMADACoherence
gapless excitation structure
spin renormalization factor

Z s =|< S = 1, q | ∑ c †↑ (k + q )c↑ (k ) − c†↓ (k + q )c↓ ( k ) | S = 0, q = 0 >|k=|< S = 1, q | ∑ S + (q ) | S = 0, q = 0 >|kZs/√N→0 for N→∞

1

Spin Renorm. / N

no fractionalization

particle-holeexcitation
Collinear
AFI0.5

Dimer
Dimmer
NMI000.20.1

0.2

1/2

1/N

M. IMADAExotic Spin Liquid at Registered Phase
specific heatmagnetic susceptibility

C2 : 7 × 7commensuratephase
2nd layer atom
1st layer atom

K. Ishida et al.,PRL 79, 3451 (1997)

T > 10µK; no order

Masutomi, Karaki andIshimoto PRL (2004)

gapless excitation
Exchange coupling ~ 10mK

M. IMADAκ−(ET)2X

X= N(CN)2Br= N(CN)2Cl= Cu(NCS)2t’
Dimerization“Approximate” single band 2D systemfrom antibonding orbitalMoriat a dimer moleculetHalf filled
anisotropic
triangular latticewith single bandtt’U

M. IMADAκ−(ET)2Cu2(CN)3

X=Cu 2 (CN) 3t '/ t ~ 1.0,
(largest among ET family)
U /t ~ 8
Shimizu, Maesato,Saito, Miyanaga, Kanoda (2003)

No signature ofmagnetic transition
down to 0.03K
No long-ranged order
Gapless excitations
M. IMADAElectronic Structure Calculationof Realistic Systems

M. IMADADownfolded HamiltonianImai, Solovyev,MI
PRL, 95 (2005) 176405

1††
H = ∑ε knc knckn +c knckn'Unn',mm' (k , k ')c k ' m ck ' m'∑
2 k ,k ',n,m,n',m'kn†

solver; PIRG-----------------------------xy=1, yz=2, zx=3------------------------------interaction(1) intra-orbital Coulomb interaction:
U(1)=2.772451, U(2)=2.583078, U(3)=2.583069(2) inter-orbital Coulomb interaction:U'(1,2)=1.346084, U'(1,3)=1.346081,U'(2,3)=1.279618(3) exchange and pair hopping:J(1,2)=0.654524, J(1,3)=0.654524, J(2,3)=0.639163

hoppingvector: 1.000000 0.000000-0.2184 0.0000 0.0000
0.0000 -0.0451 0.0000
0.0000 0.0000 -0.1940
vector: 1.000000 1.000000-0.0753 0.0000 0.0000
0.0000 0.0094 0.0029
0.0000 0.0029 0.0094
vector: 0.500000 0.500000
0.0020 0.0056 0.0056
0.0056 -0.0157 -0.0089
0.0056 -0.0089 -0.0157
vector: 2.000000 0.000000-0.0035 0.0000 0.0000
0.0000 0.0004 0.0000
0.0000 0.0000 0.0175

0.000000

0.000000

1.638780

0.000000

Renormalization factor
level-0.8468 -0.9288

degeneracy of yz,zx

renormalization factors:
ZyzM.
=0.80
Zxy=0.81
IMADA

I. Solovyev and M. ImadaSr2VO4

Similarity to high-Tc
cuprates
Layered perovskiteVO

one 3d electron：d 1 system
Nozaki et al. 1991

Srcuprates； d 9 systemone 3d hole
M. IMADA

Total energy comparison
λ:interaction parameter U⇒λU

realisticinteraction

Spin singlet ground state
PIRG result; λc~ 0.9
Close to the MI transition
Imai, Solovyev,MIPRL, 95 (2005) 176405

M. IMADAOrbital-SpinStructure

orbital stripe &plaquette spin

4×2 structure
highlycharge density
nontrivial
on the mirror planestructure

M. IMADAElectronic Structure by Our Method I
Sr2VO4

conduction

gapmagnetism

LDA0para

Hartree-Fockmetalinsulator

0.3eV

ferro

experimentslightly
insulating 〜 0-0.15eV

DFT-PIRG(close totransition）

AF?
not well known

close tonontrivial AF,
transition 〜 0-0.1eV stripe orbital(slightlyplaquette spin
insulating）

order

experimental
test?
M. IMADAElectronic structure of YVO3

OtsukaM. IMADAElectronic structure of YVO3

gap ~ 0.7 eVOtsuka

M. IMADA

Summary
Algorithm of Auxiliary-Field Path Integral forFermions
PIRG; quantum number projection
Applications
Hubbard models; Mott transitions
quantum spin liquids
Realistic models; DFT+PIRG
downfolding + low-energy solver
M. IMADA

Outlook
Convergence depends on lattice structure, interaction
and system size
Local minimum

M. IMADA