Transcript Monte Carlo

DIAGRAMMATIC MONTE CARLO:
From polarons to path-integrals to skeleton technique
N. Prokof’ev
KITPC 5/13/14
AFOSR
MURI
Advancing Research in Basic Science and Mathematics
Classical MC

Z y   d x1d x 2

d x N W x1 , x 2 ,
xN , y

the number of variables N is constant
Quantum MC (often)


Z y    d x1d x 2
n 0 

d x n Dn  ; x1 , x 2 ,
xn , y

Integration variables
term order
different terms of
of the same order
Contribution to the answer
or weight (with differential measures!)


A y    d x1d x 2
n 0


d x n Dn  ; x1 , x 2 ,

x n , y   D

Monte Carlo (Metropolis) cycle:
Accept with probability
Diagram
D
suggest a change
Same order diagrams:
Updating the diagram order:
Racc 
D ' (d x) n

O (1)
n
D (d x)
D '
D
Business as usual
D ' (d x) n  m
m


(
d
x
)

n
D
(d x)
Ooops
Balance Equation:
D 

If the desired probability density distribution
of different terms in the stochastic sum is P
then the updating process has to be stationary
with respect to P (equilibrium condition). Often
updates '
'
 ( ')R
accept 
Flux out of


updates '
P  D
'
D '  ' ()Raccept
Flux to

 ( ') Is the probability of proposing an update transforming  to  '
Detailed Balance: solve equation for each pair of updates separately
'
 '
D  ( ')R

D

(

)

R
accept
'
'
accept
Let us be more specific. Equation to solve:

Dn x1 ,


, x n (d x)n n,n m x1 ,
D 


n n  m
, x n m (d x) m Raccept
 Dn m x1,
 ( ')

m n
, x nm (d x) n m n m,n Racncept
D 
  ' ( )
new variables x n 1 ,
, xnm
are proposed from the
normalized probability distribution
Solution:
R
n  m
Rnaccept

n  m n
Raccept
Dn m ( x1 ,
Dn ( x1 ,
n  m , n
, xn m )

, xn ) n,n m ( x1 , , xn m )
All differential measures are gone!
Efficiency rules:
- try to keep R ~ 1
- simple analytic function
n,n  m ( xn 1 ,
, xn  m )
ENTER
Standard Monte Carlo setup:
- configuration space
(depends on the model
and it’s representation)
- each cnf. has a weight factor
 Ecnf / T
cnf
e
W
A W

W
cnf
- quantity of interest
Acnf
A
cnf
cnf
cnf
cnf
Statistics:

e
 E state / T
MC
O state
{ states }
Monte Carlo

Ostate 
{ states }
states generated from
probability distribution
Anything:

Monte Carlo
i arg[ F ( )]
e
O ( )

{ }
states generated from
probability distribution
Connected Feynman
Anything = diagrams, e.g. for the
proper self-energy

MC
F ( )O ( )
{  any set
of variables}
diagram order
topology
internal variables
e  E state / T
Monte Carlo
| F ( ) |
Answer to S. Weinberg’s
question
MC
   sign ( F (v ))

Configuration space = (diagram order, topology and types of lines, internal variables)
Diagram order
{qi , i , pi }
Diagram topology
This is NOT: write/enumerate diagram after diagram,
compute its value, and then sum
Polaron problem:
H  H particle  H environment  H coupling 
quasiparticle
E ( p  0), m* , G( p, t ), ...
Electrons in semiconducting crystals (electron-phonon polarons)
H    ( p)a p a p 
electron
p
e
e


(
p
)(
b

q bq  1/ 2) 
phonons
q
 V a
q
pq

p q
a p bq  h.c.
el.-ph.
interaction
H    ( p)a p a p    ( p)(bqbq  1/ 2)   Vq a p q a p bq  h.c.
p
q
pq
electron
phonons
el.-ph. interaction
q
q
p
pq
p
Green function:
G( p, )  a p (0)a  p ()  a p e - H a  p e  H

p

0
+
+
Sum of all Feynman diagrams
Positive definite series in the
+
( p, )
representation
+
+
+ …
q2
q1
G( p, )  
Feynman
digrams
q2
p  q2
p  q1
p
0
q3
1
p
1 '
4 '

Graph-to-math correspondence:
 


G p,    d x1d x 2
d x n Dn  ; x1 , x 2 ,
n 0 

x n , p, where xi =(qi , i , i ')
is a product of
qi
qi
1
Vqi
Positive definite series in the
e
pi
1 '
  ( pi )( i '  i )
( p, ) representation
1
1 '
e  ( q )( i ' i )
0
Diagrams for:
bq1 (0)bq2 (0)a p  q1  q2 (0)a p  q1  q2 ( )bq1 ( ) bq2 ( )
there are also diagrams for optical conductivity, etc.
Doing MC in the Feynman diagram
configuration space is an endless fun!
The simplest algorithm has three updates:
Insert/Delete pair (increasing/decreasing the diagram order)
q2
q1
p1  p
p3  p
p2
1 '
1
D
p1
p2 '
2
q1
2 '
1
p5  p
p4 '  p2
p3 '
1 '
D '
D ' / D  | Vq2 |2 e ( q2 )( 2 ' 2 ) e( ( p2 ') ( p2 ))(1  2 ) e( ( p3 ') ( p2 ))( 2 '1 )
n1,n
D '
D '
1
R 


D n,n 1 ( x1 , , xn 1 ) D (n  1)n,n1 ( x1 ,
In Delete select the phonon line to be deleted at random
, xn1 )
The optimal choice of
Frohlich polaron:
2
 n ,n 1 ( x1 ,
, xn 1 )
  p2 / 2m, q  0 , Vq   / q
q 0 ( 2 ' 2 )
D ' / D  2 e
e
q
[( p2 ')2  p22 ](1  2 ) [( p3 ')2  p22 ]( 2 ' 1 )

2m
1. Select  2 anywhere on the interval
2. Select  2 ' anywhere to the left of
3. Select
depends on the model
2
(0, ) from uniform prob. density
 ( '  )
from prob. density e 0 2 2
(if  2 '   reject the update)
q2 from Gaussian prob. density e
 n ,n 1 ( 2 , 2 ', q2 ) e
dqd d d 2
 ( q22 / 2 m )( 2 '  2 )
, i.e.
0 ( 2 '  2 )  ( q22 /2 m )( 2 '  2 )
e
New
:
p
 last 
 '
Standard “heat bath” probability density
Always accepted,
~ e  ( p )( ' last )
R 1
Normalization:
histogram
hbin 

f ( x)dx
bin
special “bin” where

f ( x)dx  I
bin
is known exactly
I
x
Normalized histogram
hbin
I
hsbin
I
x
Normalization using “desined bin”:
hbin 

f ( x)dx
bin
I0
x
hbin
I0
h0
I0
x
This is it! Collect statistics for
Analyze it.
G ( p , )
or some other corr. function.
Analysing data:
G( p,   )  Z peE ( p)
Z p | C p |
2
dispersion relation
probability of getting
a bare electron
(Lehman expansion)
E p  C p a †p 0   C p ,q bq† a †p q 0   C p ,q1 ,q2 bq†1 bq†2 a †p q1 q2 0  ...
q
ln G p ( )
q1q2
Z p(2)   | C p,q1 ,q2 |2
q1 , q2
1
Zp
 Ep
01

probability of getting
two phonons in the
polaron cloud
 0
High-energy
Standard model
physics
High-Tc
Hubbard
superconductors
model
Quantum chemistry
Coulomb&gas
band structure
Quantum
Heisenberg
magnetism
model
Periodic Anderson
Heavy fermion
& Kondo
materials
lattice models
…
- Introduced in mid 1960s or earlier
- Still not solved (just a reminder, today is 01/13/2012)
- Admit description in terms of Feynman diagrams
Feynman Diagrams & Physics of strongly correlated many-body systems
In the absence of small parameters, are they useful in higher orders?
series areskeleton
the devil'sgraphs!
invention...”
for regularized
N.Abel,Yes,
1828:with sign-blessing“Divergent
And if they are, how to handle millions and billions of skeleton graphs?
them
withToday,
Diagrammatic
StevenSample
Weinberg,
Physics
Aug. 2011 Monte
: “Also,Carlo
it was techniques
easy to imagine
any number
of quantum field
(teach computers rules of quantum
field theory)
theories of strong interactions
but what could anyone do with
Unbiased
From current strong-coupling
them?” solutions beased
on millions of graphs with
theories based on one lowest
extrapolation to the infinite
order skeleton graph (MF,
diagram order
RPA, GW, SCBA, GG0, GG, …
Skeleton diagrams up to high-order: do they make sense for g  1 ?
NO
YES
g even if they
Series diverge for large
converge for g  1
Divergent series outside of finite
convergence radius can be re-summed.
n
Dyson “Expansion in g is asymptotic if for
some g the system becomes pathological”:
Continuous space bosons and fermions
collapse to infinite density for U   U
Math. Statement: number of skeleton graphs
 2n n3/2 n !
asymptotic series
with zero conv. radius
n
Skeleton series are not based on g
Many systems remain well-defined for g
Lattice fermions, quantum magnetism,
resonant fermions … _+ regularization
Number of graphs is  2 n n !
but because they alternate in sign they
may all cancel each other to near zero !
n
3/2
“Sign blessing”
AN
Asymptotic series
for g  1 with zero
convergence radius
are useless!
Start computing high-order diagrams!
1/ N
Re-summation of divergent series with finite convergence radius.

Example:
A   cn  3  9 / 2  9  81/ 4  ... 
бред какой то
n 0
Define a function
such that:
f n, N
f n, N
Nb
 1 for n  N 1
f n, N  e
 n2 / N
f n, N  e n ln( n )
f n , N  0 for n  N
Na
n

Construct sums
AN   cn f n , N
n 0
and extrapolate
lim AN
N 
to get
AN
ln 4
1/ N
A
(Gauss)
(Lindeloef)
Conventional Sign-problem vs Sign-blessing
Computational complexity is exponential in system volume
Sign-problem:
(diagrams for
Z)
tCPU  exp{# Ld  }
and error bars explode before a reliable
exptrapolation to L   can be made
Feynamn diagrams: No
(for
ln Z )
limit to take, selfconsistent formulation, admit
analytic results and partial resummations.
Sign-blessing:
(diagrams for
L 
ln Z
Number of diagram of order
)
n
is factorial
 n !2n n3/2 thus the
only hope for good series convergence properties is sign alternation of diagrams leading to their cancellation. Still,
tCPU  n !2n n3/2
i.e. Smaller and smaller error bars are likely
to come at exponential price (unless convergence is exponential).
Diagrammatic Monte Carlo in the generic many-body setup
Feynman diagrams for free energy density
G ( p , )
x
U (q )
G ( p , )
x
x
x
G
U


(0)

G
U
+
-
G(0)
U


G
U
Bold (self-consistent) Diagrammatic Monte Carlo
Diagrammatic technique for ln(Z ) diagrams: admits partial summation and
self-consistent formulation
No need to compute all diagrams for G and
G ( p, )

Dyson Equation:
+
Gp
( ,)
0
U:
(
p,1   2 ) +

Screening:
U
+
+ ...

+
U
Calculate irreducible diagrams for  ,  , … to get G ,
U,
…. from Dyson equations


+
+
+
. . +. . +.
+
. +.
.+
In terms of “exact” propagators
Dyson Equation:

+

Screening:

+

G(0)
More tools: Build diagrams using ladders:
(contact potential)

(0)

+
 [
-
+
U
G(0)
+
] +
+
In terms of “exact” propagators
Dyson Equations:

+


+

Fully dressed skeleton graphs (Heidin):

+




+



Irreducible 3-point vertex:
3 

3
.

.

1
all accounted for already!
 ...
Unpolarized system at unitarity:BCS-BEC crossover
Unitary gas:
kF aS  
when
kF
and
F
are the only length/energy scales
Answering Weinberg’s question:
Equation of State for ultracold fermions & neutron matter at unitarity
MIT group: Mrtin Zwierlein, Mark Ku, Ariel Sommer,
Lawrence Cheuk, Andre Schirotzek
Uncertainty due to location of
the as   resonance B  834 1.5G
BDMC results
Kris Van Houcke, Felix Werner,
Evgeny Kozik, Boris Svistunov, NP
virial expansion (3d order)
Ideal Fermi gas
QMC for connected Feynman diagrams NOT particles!
Sign blessing
Sign problem
H  H 0  H1   U ij ni n j   i ni   t (ni , n j )b j bi
ij
ij 
i
Lattice path-integrals for bosons and
spins are “diagrams” of closed loops!

Z  Tre
H
TreH0
 Tre
H 0
e

 H1 (  )d 
0

 




1   H1 ()d     H1 ()H1 ()d d  '  ...


0 
 0


imaginary time
t (1, 2)
'
+
0
i
j

i
j
+
ni  0,1, 2,0
Diagrams for
Diagrams for
Z = Tr e-  H


imaginary time
imaginary time
GIM =Tr T bM† ( M ) bI ( I ) e-  H
M
I
0
0
lattice site
lattice site
The rest is conventional worm algorithm in continuous time
Simulating Bose-Hubbard model “as is” and comparing to experiments:
(in this example, N  300000)
Nature Physics, 6, 998–1004, (2010)
It is realistic to do about 2,000,000 or more particles at temperatures
relevant for the experiment.
Phys. Rev. Lett. 103, 085701 (2009)
2


 P  /   m( Ri 1  Ri )

Z   dR1...dRP exp    
U ( R)  
2


 i 1 

Path-integrals in continuous space
are “diagrams” of closed loops too!
P

ri ,1 
2
1
P
ri ,2 
Ri  (ri ,1 ,ri ,2 ,...,ri , N )
Diagrams for the attractive tail in U ( r ) :
If 
N
U ( r
j i
j
 ri )  (| rj  r | rc )  1 and N>>1 all the effort is for something small !
e
V ( rij )
 1  (e
V ( rij )
 1)  1  pij
M asha
statistical interpretation
Ira
ignore V ( rij )
p ij
i
j
: stat. weight 1
Account for V ( rij ) : stat. weight p
Faster than conventional scheme for N>30, scalable (size independent) updates
with exact account of interactions between all particles
3D @ s.v.p.
64
experiment
2048
TCAziz  2.193vs TCexp  2.177
LS / T  W 2 / 3
Other applications: Continuous-time QMC solves (impurity solvers)
are standard DMC schemes
Fermions with contact interaction
1
3
G11 G12 ... G1 p
2
det Gij 
9
  {n; r1 ,1;r2 , 2 ;...;
  rn , n }
G21 G22 ... G2 p
... ... ... ...
G p1 G p 2 ... G pp


Z    ... ( d rd ) p ( U ) p det G ( x i , x j ) d et G ( x i , x j )
Rubtsov (2003)
p 0
Most efficient solvers for DMFT and DCA are based on this approach
0.50
M. Holland et al., PRL 87, 120406 (2001)
Experiment
Analytical results
Numerical results
This work
0.30
J. Kinhast et.al. Science 307 1296 (2005)
0.25
A. Bulgac et al., cond-mat/0505374
Tc/EF
0.20
P. Nozieres and S. Schmitt-Rink, J. Low. Temp. Phys 59, 195 (1985)
0.15
0.10
Maier et al.
A. Sewer et al.,
PRB 66, 140504 (2002)
0.05
X.-J. Liu and H. Hu, cond-mat/0505572
M. Wingate, cond-mat/0502372
0.00
0.0
0.1
0.2
0.3
0.4
0.5

L3  20 3 , N el  250, F
0.1
1/3
0.6
0.7
0.8
0.9
1.0
Critical point from pair distribution function
Criticality:
(k  0)   2  (T  TC )  (2 )
(k ,   0)
(  )C
from zero of
  0.038  0.6717
A  (k  0) 1/(2 )  T  TC
TC
  2.25
1

 F 2