Computational methods in Coupled Electron-Ion Monte Carlo

0.22 MB·

by  Carlo Pierleoni

5
0
0
2
Computational methods in Coupled Ele
tron-Ion

n
a
J
Monte Carlo
4

] a ∗ b
h
Carlo Pierleoni and David M. Ceperley
p
-
p
m 9th February 2008
o
c
. a
s
c Department of Physi
s, University of L'Aquila, Via Vetoio, I-67010 L'Aquila,
i
s Italy
y
b
h
Dept. of Physi
s and NCSA, University of Illinois at Urbana-Champaign, Ur-
p
[
bana, IL 61801, USA

1
v
Abstra
t
3
1
0 In the last few years we have been developing a Monte Carlo simula-
1
0 tion method to
ope with systems of many ele
trons and ions in the Born-
5
Oppenheimer (BO) approximation, the Coupled Ele
tron-Ion Monte Carlo
0
/
s Method (CEIMC). Ele
troni
properties in CEIMC are
omputed by Quan-
c
i tum Monte Carlo (QMC) rather than by Density Fun
tional Theory (DFT)
s
y
based te
hniques. CEIMC
an, in prin
iple, over
ome some of the limita-
h
p
tions of the present DFT based ab initio dynami
al methods. Appli
ation of
:
v
thenewmethodtohighpressuremetalli
hydrogenhasre
entlyappeared. In
i
X
thispaperwepresentanewsamplingalgorithmthatwehavedevelopedinthe
r
a
framework of the Reptation QuantumMonte Carlo (RQMC) method
hosen
to sample the ele
troni
degrees of freedom, thereby improving its e(cid:30)
ien
y.
Moreover,weshowherethat,atleastforthe
aseofmetalli
hydrogen,varia-
tionalestimatesof theele
troni
energies leadto ana

uratesampling ofthe
protondegrees of freedom.
∗
Corresponding author: Carlo Pierleoni, Physi
s Dept. University of L'Aquila, Via Vetoio,
67010Coppito,L'Aquila(Italy),fax:[phone];email:
arlo.pierleoniaquila.infn.it
1
1 Introdu
tion
Modern theoreti
al methods in
ondensed matter physi
s and
hemistry rely heav-
ily on numeri
al simulations. The problem of solving the S
hroedinger equation
for many-body systems is too di(cid:30)
ult to be addressed dire
tly, even within the
simpli(cid:28)
ationprovidedbytheBorn-Oppenheimerapproximation. Inthemostpop-
ular pra
ti
al approa
hes (Hartree-Fo
k (HF) and the Density Fun
tional Theory
(DFT)basedmethods[1℄)theoriginalproblemisrepla
edbytheproblemofsolving
the time independent S
hroedinger equation for a single ele
tron in the (cid:28)eld of the
nu
leiandthemean (cid:28)eldgeneratedbytheotherele
trons. DFT is, in prin
iple,an
exa
ttheory but the energyfun
tional mustbe treated approximatelyforpra
ti
al
purposes. In the simplest Lo
al Density Approximation (LDA), this exa
t theory
be
omes a self-
onsistent mean (cid:28)eld theory. Extensions of LDA, su
h as Gener-
alized Gradient Approximation (GGA) provide more a

urate results but remain
essentially at the level of an e(cid:27)e
tive mean (cid:28)eld treatment. Despite the mean (cid:28)eld

hara
ter, DFT s
hemes have proved to provide quite a

urate results for many
di(cid:27)erent systems[1℄
In 1985, Car and Parrinello introdu
ed an e(cid:30)
ient method to
ouple standard
Mole
ular Dynami
s for
lassi
al nu
lei with the ele
troni
stru
ture
al
ulation at
the level of LDA done (cid:16)on the (cid:29)y(cid:17) to extra
t the nu
lear for
es[2℄. Be
ause the
method allowed study of the statisti
al me
hani
s of
lassi
al nu
lei with many
body ele
troni
intera
tions, it opened the way for the use of simulation methods
for realisti
systemswith an a

ura
ywell beyond the limits ofe(cid:27)e
tive for
e (cid:28)elds
available. Inthelasttwentyyears,thenumberofappli
ationsoftheCar-Parrinello
ab-initio mole
ular dynami
s has ranged from simple
ovalent bonded solids, to
highpressurephysi
s,materials
ien
eandbiologi
alsystems. Therehavealsobeen
extensionsoftheoriginalalgorithmtosimulatesystemsat
onstanttemperatureand

onstant pressure[3℄, (cid:28)nite temperature e(cid:27)e
ts for the ele
trons [4℄, and quantum
nu
lei [5℄.
Despite re
ent progress, DFT su(cid:27)ers from well-known limitations, for example,
ex
ited state properties su
h as opti
al gap and spe
tra are less reliable. DFT
showsseriousde(cid:28)
ien
iesin des
ribingvan derWaalsintera
tions, non-equilibrium
geometries su
h as rea
tion barriers, systems with transition metals and/or
luster
isomerswith
ompetingbonding patterns[1,6℄. Asa
onsequen
e,
urrentab-initio
predi
tionsofmetallizationtransitionathighpressures,orevenpredi
tionofphase
transitionsareoftenonlyqualitative. Hydrogenisanextreme
ase[7,8,9℄ but even
β
in sili
on the diamond/ -tin transition pressure and the melting temperature are
seriously underestimated[10℄.
Another route to the ground state properties of a system of many ele
trons
in presen
e of nu
lei is the Quantum Monte Carlo method[11, 6℄. In its simplest
form, an analyti
many ele
tron wave fun
tion is
hosen on the basis of the vari-
ational prin
iple (Variational Monte Carlo, VMC) and the quantum averages are
obtained by a Metropolis Monte Carlo simulation of the ele
troni

oordinates. A
morea

uraterepresentationofthe groundstatewavefun
tion
anbeobtained by
exp{−β H} H
e
proje
ting the variational wave fun
tion with the operator where
β
e
is the many-body hamiltonian, and is the proje
tion time. Provided that the
variational wave fun
tion is not orthogonal to the ground state wave fun
tion, the
proje
ted fun
tion tends exponentially fast to the ground state wave fun
tion as
β → ∞
e
. Sin
e matrix elements of the above proje
tion operator at large values
β P
e
of are unknown for non trivial systems, a Trotter breakup in many ( ) small
τ =β /P)
e e
imaginary time intervals ( must be employed. In the
on(cid:28)guration rep-
3N
resentation, ea
h proje
tion
orresponds to a -dimensional integral whi
h
an
be performed by Metropolis Monte Carlo method provided that the propagator in
imaginarytime
anbe
hosenrealand
anbe interpreted asaprobabilitydistribu-
tion. This is the essen
e of the Di(cid:27)usion Monte Carlo method (DMC) whi
h is an
(cid:16)exa
t(cid:17) methodforsystemsofbosonsorboltzmannons. Thismeansthatallsystem-
ati
errors in a simulation are under
ontrol in the sense that they
an be redu
ed
asmu
hasdesired. Sin
e ele
tronarefermions,the aboves
heme fails be
ausethe
imaginary time propagator must be
ompletely antisymmetri
under ex
hange of
two ele
trons and therefore
annot be
hosen stri
tly non-negative everywhere in

on(cid:28)gurational spa
e. This is the origin of the infamous (cid:16)fermion sign problem(cid:17).
In order to avoid the sign problem the (cid:16)(cid:28)xed node approximation(cid:17) has been pro-
posed and used routinely to perform fermion simulations[6℄. The energy
al
ulated
withthis approximationisvariationalwithrespe
ttothepositionofthenodalsur-
fa
es of the trial wave fun
tion. Over the years, the level of a

ura
y of the (cid:28)xed
3He
nodeapproximationforsimplehomogeneoussystems,su
has andtheele
tron
gas,hasbeensystemati
allyimprovedbyintrodu
ingmoresophisti
atednodalsur-
fa
es (ba
k(cid:29)ow orbitals)[12, 13℄. In more
omplex, inhomogeneous situations su
h
as atoms, mole
ules and extended systems of ele
trons and nu
lei, progress have
been somewhat slower. Nonetheless, in most
ases, (cid:28)xed-node QMC methods have
provedtobemorea

uratethanmean(cid:28)eldmethods(HFandDFT)[6℄. Computing
ioni
for
eswith QMCto repla
ethe DFT for
esin theab-initioMD, ismoredi(cid:30)-

ult and a general and e(cid:30)
ient algorithm is still missing. Moreover, the
omputer
time requiredforaQMCestimateofthe ele
troni
energyis, in general,morethan
for a
orresponding DFT-LDA
al
ulation. These problems have seriously limited
the development of an ab-initio simulation method based on the QMC solution of
the ele
troni
problem (cid:16)on the (cid:29)y(cid:17).
In re
ent years, we have developed a di(cid:27)erent strategy based entirely on the
Monte Carlo method both for solving the ele
troni
problem and for sampling the
ioni

on(cid:28)gurationspa
e[14,15℄. Thenewmethod,
alledtheCoupledEle
tron-Ion
Monte Carlo method (CEIMC) has been applied so far to high pressure metalli
hydrogen where it has found quite di(cid:27)erent e(cid:27)e
ts of temperature than CPMD
basedontheLDAfor
es[16℄. Ourpresentinterpretationofthedisagreementisthat
LDA providesaBorn-Oppenheimersurfa
e quitesmootherthanthe morea

urate
T >0
QMC one and this strongly a(cid:27)e
ts the stru
ture of the protoni
system at .
Thepaperisorganizedasfollows. Thefollowingse
tion2isdevotedtoanoutline
of the CEIMC method. We will not go into all details sin
e two long arti
les have
appearedon generalaspe
ts and earlyimplementationsofthe method[14,15℄. One
of the new aspe
ts that we havere
ently implemented in CEIMC, not des
ribed in
thosereferen
es,istheReptationQuantumMonteCarloproje
tionoftheele
troni
variationalwavefun
tion[17℄. Sointhesubse
tion2.1wereviewtheRQMCmethod
andinthefollowingsubse
tion2.2wefo
usonthesamplingalgorithm,weintrodu
e
our new s
heme to improve e(cid:30)
ien
y and reliability of RQMC and we provide an
analyti
al proof. In Se
tion3 we report numeri
al results on the
onvergen
eof the
new s
heme with the proje
tion time and with the Trotter time step. Finally, in
se
tion 4, we
on
lude.
2 The Coupled Ele
tron-Ion Monte Carlo method
CEIMC method is based on the Born-Oppenheimer (BO) separation between the
slow nu
lei and the fast ele
trons. This is in
ontrast with other Quantum Monte
Carlo methods, Di(cid:27)usion Monte Carlo (DMC)[11, 6℄ or (cid:28)nite temperature Path
IntegralMonte Carlo(PIMC)[18,19℄ methods where ele
tronsand ions aretreated
on the same footing. As usual, the BO approximation allows to over
ome the
limitations of the other QMC methods, while introdu
ing an often negligible error.
InCEIMC,the
on(cid:28)gurationalspa
eoftheprotondegreesoffreedomatinverse
β = (k T)−1
B
temperature is sampled with a Metropolis algorithm in whi
h the
S S′
di(cid:27)eren
e between the BO energy of a proton state and of a trial state is

omputed by an ele
troni
ground state QMC
al
ulation. The QMC estimate of
∆=[E(S′)−E(S)]
the energydi(cid:27)eren
e hasstatisti
alnoisewhi
hwouldbias the
standard Metropolis algorithm. Unbiased sampling of the proton
on(cid:28)gurations is
∆
a
hieved by the penalty method[20℄ whi
h repla
es the energy di(cid:27)eren
e in the
∆+(βσ∆)2/2 σ∆2
a

eptan
eformulaby ,where isthevarian
eoftheenergydi(cid:27)er-
σ2 > 0
∆
en
e. Sin
e , the noise always
auses extra reje
tions but this
ompensates
for (cid:17)uphill(cid:17) moves a

epted be
ause of a favorable energy (cid:29)u
tuation.
Severalmethodsfor
omputingenergydi(cid:27)eren
esinQMCareavailable.[14,15℄.
Asimpleand e(cid:30)
ientmethodisto sampletheele
troni
degreesoffreedom from a
distributionfun
tionwhi
histhesumoftheele
troni
distributionfun
tionsforthe
S S′
and states(e. g. thesumofthethesquareofthetrialwavefun
tionsinVMC).
Averages of operators involving ele
troni
degrees of freedom and a single proton
S

on(cid:28)guration, say , (for instan
e total energy, varian
e, et
.) are then
omputed
by
orrelatedsampling[11,14,15℄. Forthetypi
alsizeoftheprotonmoves(between
0.01 0.5
Å and Å for
lassi
al protons depending on density and temperature) and
thetypi
alsystemsize(upto54protons)wehaveinvestigated,thismethodismu
h
moree(cid:30)
ientthanperformingtwoindependentele
troni

al
ulationsforthestate
S S′
and .
In the ground state QMC methods, an ele
troni
trial wave fun
tion must be

hosena

ordingtothephysi
softhesystembeingstudied. Forthemetalli
phase
of hydrogen, we have re
ently developed analyti
fun
tions whi
h in
lude ba
k(cid:29)ow
and three-body
orrelations[21℄. These wave fun
tions are parti
ularly appropriate
to CEIMC sin
e they have a

urate energies already at the variational level, they
havenoadjustableparametersrequiringoptimization,andtheir
omputational
ost
is mu
h less than using orbitals expanded in a plane wave basis typi
ally used in
QMC
al
ulations[22℄.
Inmetalli
systems,(cid:28)nitesizee(cid:27)e
tsarelargeandmustbesuitablytreated. The

ommon pro
edure is to repeat the
al
ulation for systems of in
reasing size and
extrapolatetothethermodynami
limitbutthisisimpra
ti
alwithinCEIMC,sin
e
itwouldhavetobeperformedforanyproposedprotoni
stepbeforeitsa

eptan
e.
A mu
h better strategy is to use Twist Averaged Boundary Condition (TABC)
1/N
[23,15℄whi
hredu
esthe(cid:28)nitesizeerrorintheenergytothe
lassi
al behavior.
It
onsistsinaveragingtheenergyoverthephasethatthemanybodywavefun
tion

an pi
k if a single ele
tron wraps around the super-
ell. This is equivalent to
Brillouin zone sampling in the single ele
tron approximation. Within CEIMC it
does not
ause a large in
rease in required CPU time/step.
Finally a re
ent improvement of the method is the introdu
tion of quantum
e(cid:27)e
ts for the protons, quite important in high pressure hydrogen. This is done by
developing the thermal density matrix of protoni
degrees of freedom on the BO
surfa
e in Feynman Path Integrals[18, 16℄. A similar te
hnique in the
ontext of
Car-Parrinello method has appeared[5℄. We are not going to dis
uss the last two
aspe
tsoftheCEIMC.WhileTABCimplementationinCEIMChasbeendes
ribed
in ref.[15℄, our implementation of PIMC for proton degrees of freedom in CEIMC
will be the subje
t of a future publi
ation.
2.1 Reptation Quantum Monte Carlo Method
To go beyond VMC ele
troni
energies, we implemented a Reptation Quantum
Monte Carlo algorithm (RQMC)[17℄, rather than Di(cid:27)usion Monte Carlo algorithm
(DMC). The implementation of the energy di(cid:27)eren
e method is more straightfor-
ward in RQMC, nor are averages of observables whi
h do not
ommute with the
hamiltonian biased.
In RQMC the ground state wave fun
tion is obtained by
onstru
ting an imag-
|Ψ0i
inary time path integral for the ele
troni
degrees of freedom. If is the trial
state,thetrialstateproje
tedina(cid:16)time(cid:17) βe/2,|Ψβe/2i=e−βeH/2|Ψ0i,will
onverge
β
e
to the ground state for large . Let us de(cid:28)ne the (cid:16)partition(cid:17) fun
tion
Zβe =hΨ0|e−βeH|Ψ0i. (1)
The energy is then de(cid:28)ned as
d 1
E(βe)=−dβ lnZβe = Z hΨ0|e−βeHH|Ψ0i=hEL(R)i (2)
e βe
EL(R) = ℜ(Ψ−01(R)HΨ0(R))
where the averages of the lo
al energy, , are with
respe
t to the path average. In pra
ti
e, the energy is
omputed as the average
ℜ
of the lo
al energy at the two ends of the path. Here indi
ates the real part in
E(β )
e
the
ase that the trial fun
tion or Hamiltonian is
omplex. The energy is an
β
e
upper bound to the (cid:28)xed node energy for ea
h value of , it
onverges to this at
β β
e e
large , and its derivative is stri
tly negative. This latter quantity is, in fa
t,
minus the varian
e of the total energy
dE(β )
σ2(β )=− e =hE (0)E (β )i−hE(β )i2.
e L L e e
dβ (3)
e
β
e
The varian
e tends to zero for large enough values of providing a useful signal
for the
onvergen
e of the energy to the ground state. This is the zero varian
e
theorem in RQMC. On the other hand, the varian
e of the lo
al energy
omputed
σ2(β )
e
at either end of the path is the mixed estimator of DMC for . In pra
ti
al
β
e
implementations, it is desirable to keep as small as possible to maximize the
e(cid:30)
ien
y of the energy di(cid:27)eren
e method.
To
ompute the needed density matrix elements, we divide the proje
tion time
β P τ = β /P
e e e
into time sli
es and make a semi-
lassi
al approximation for
exp(−τeH) R={r1,...,rN}
. Our notation for a single ele
troni

on(cid:28)guration is ,
s = {R0,R1,...,RP}
while for the entire path is . The probability distribution for
a path is P−1
Π(s)=exp −U(R0)−U(Rp)− Ls(Ri+1,Ri)
( ) (4)
i=1
X
U(R) = ℜ[lnΨ0(R)] Ls(R,R′)
where and is the symmetrized link a
tion for our
approximationoftheshorttimepropagator. Wehaveusedtheimportan
esampling
Green's fun
tion of the DMC propagator
R|e−τeH|R′ = Ψ0(R) exp −τ E (R)− [R′−R−2λτeF(R)]2 .
(cid:12)Ψ0(R′)(cid:12) (cid:20) e L 4λτe (cid:21) (5)
(cid:10) (cid:11) (cid:12) (cid:12)
(cid:12) (cid:12)
(cid:12) (cid:12) 2
F(R)=∇U(R) λ=h¯ /2m
e
wherethefor
eis and ,whi
hprovidesthesymmetrized
link a
tion
τ
L (R,R′) = e E (R)+E (R′)+λ F2(R)+F2(R′)
s L L
2
(R(cid:2)−R′)2 (R−R′)·(F(cid:0)(R)−F(R′)) (cid:1)(cid:3)
+ +
4λτ 2 (6)
e
An alternative form for the link a
tion
ould be obtained through the pair a
tion
developed in (cid:28)nite temperature Path Integral MC[18℄. However, we have not im-
plemented this form and do not have a
omparison of its e(cid:30)
ien
y.
Inordertoimposethe(cid:28)xedphase
onstraintontheproje
tedwavefun
tion,we
LFP(R,R′)=λτ 1dη|∇φ(X(η))|2
mustaddtothelinka
tionatermoftheform s e 0
φ(X) R X
where isthe phase of the trial wavefun
tion at ele
troni
position and the
X(η) X(0)= R, X(1)=
integral is taken over all paths with boundary
onditions
R′
. We have taken an end-point approximation for this term ex
ept for real wave
fun
tions in whi
h
ase (cid:28)xed-node boundary
onditions were used.
Note that in the expressions above, the dependen
e on the nu
lear degrees of
freedomwasnotshowneventhoughallquantitiesdependonthem. Theprobability
Π(s,S) S
distribution of an ele
troni
path will be where indi
ates the position of
all nu
lei. Be
ause we have an expli
it distribution of the ele
troni
paths, it is
straightforwardtoapplytheimportan
esamplings
hemefortheenergydi(cid:27)eren
es
[Π(s,S)+Π(s,S′)] S S′
bysamplingtheprobabilitydistribution where and arethe
Π(s,S) ∝

urrent and the trial protoni
state, respe
tively. Note also, for VMC
2
|Ψ0(s,S)|
be
omes the squareofmodulus ofthetrial wavefun
tion(noproje
tion
R0 =RP
and ).
2.2 The (cid:16)boun
e(cid:17) algorithm
In the original work on RQMC[17℄, the ele
troni
path spa
e was sampled by a
reptation algorithm, an algorithm introdu
ed to sample the
on(cid:28)gurational spa
e
of linear polymer
hains. The slithering snake or reptation method seems to have
s
originated by Kron[24℄ and by Wall and Mandel[25℄. Given a path
on(cid:28)guration ,
R0 RP
amoveisdoneintwostages. Firstoneofthetwoends(either or )issampled
R
g
with probability 1/2 to be the growth end . Then a new point near the growth
R +2λτ F(R )
g e g
end is sampled from a Gaussian distribution with
enter at . In
orderto keep the number oflinks on the path length
onstant, the oldtail position
isdis
ardedinthetrialmove. Themoveisa

eptedorreje
tedwiththeMetropolis
formula based on the probability of a reversemove. For use in the following, let us
d d=+1 R = R d=−1
g P
de(cid:28)ne the dire
tion variable as for a head move ( ), and
Rg =R0 d
foratailmove( ). Instandardreptation,thedire
tion is
hosenrandomly
at ea
h attempted step.
P(s→s′)
In the standardreptation algorithm, the transition probability is the
T (s→s′) a (s→
d d
produ
tofanattemptprobability andana

eptan
eprobability
s′)
. Note that the path distribution given in Eq.(4), be
ause of the symmetrized
d
link a
tion does not depend on the dire
tion in whi
h it was
onstru
ted. In the
Metropolis algorithm, the a

eptan
e probability for the attempted move is
Π(s′)T (s′ →s)
a (s→s′)=min 1, −d
d Π(s)T (s→s′) (7)
(cid:20) d (cid:21)
P (s→s′)
d
whi
h ensures that the transition probability satis(cid:28)es detailed balan
e
Π(s)P (s→s′)=Π(s′)P (s′ →s)
d −d
(8)
The auto
orrelation time of this algorithm in Monte Carlo steps, that is the
[(β /τ )2/A]
e e
numberofMCstepsbetweentwoun
orrelated
on(cid:28)gurations,s
alesas ,
A β
e
where is the a

eptan
e rate, an unfavorable s
aling for large . Moreover the
o

asionalappearan
eofpersistent
on(cid:28)gurationsboun
ingba
kandforthwithout
reallysamplingthe
on(cid:28)gurationspa
ehasbeenpreviouslyobserved[26℄. Theseare
two very unfavorable features, parti
ularly in the present
ontext, where we need
to perform many di(cid:27)erent ele
troni

al
ulations (at least one per protoni
move).
There is a premium for a reliable, e(cid:30)
ient and robust algorithm.
We have found that a minimal modi(cid:28)
ation of the reptation algorithm solves
both of these problems. The idea is to
hose randomly the growth dire
tion at the
beginning of the Markov
hain, and reverse the dire
tion upon reje
tion only, the
(cid:16)boun
e(cid:17) algorithm. Asfarasweareaware,(cid:17)boun
e(cid:17) dynami
shasnotbeen previ-
ously investigated for RQMC, though Wall and Mandel[25℄ mentioned it without a
detailedproofandsubsequentpolymersimulationsdidnotuseboun
e,perhapsbe-

ausethe a

eptan
e ratioin the polymer systems is mu
h smallerthan in RQMC.
There is a related algorithmfor dire
ted loop algorithm on the latti
e and for sim-
ulations of trapped di(cid:27)usion[28, 27℄.
What follows is the proof that the boun
e algorithm samples the
orre
t prob-
Π(s) d
ability distribution . The variable is no longer randomly sampled, but,
as before, the appropriate move is sampled from the same Gaussian distribution
T (s→s′)
d
and a

epteda

ordingtothe Eq. (7). Tobeable tousethete
hniques
d
of Markov
hains, we need to enlarge the state spa
e with the dire
tion variable .
{s,d}
In the enlarged
on(cid:28)guration spa
e , let us de(cid:28)ne the transition probability
P(s,d → s′,d′)
of the Markov
hain. The algorithm is a Markov pro
ess in the
extended path spa
e, and it is ergodi
as DMC method, hen
e, it must
onvergeto
Υ(s,d)
a unique stationary state, satisfying the eigenvalue equation:
Υ(s,d)P(s,d→s′,d′)=Υ(s′,d′).
(9)
s,d
X
Π(s)
We show that our desired probability is solution of this equation. Within the
P(s,d→s′,d′)6=0 d=d′
imposed rule not all transitions are allowed, but for and
s 6= s′ d′ = −d s = s′
(a

epted move), or and (reje
ted move) only. Without loss
d′ = +1 ±1
of generality let us assume sin
e we have symmetry between . Eq. (9)
Υ(s,d) Π(s)
with repla
ed by is
Π(s′)P(s′,−1→s′,1)+ Π(s)P(s,1→s′,1)=Π(s′).
s6=s′
X
Π(s)P(s,1→s′,1)=Π(s′)P(s′,−1→
Be
auseof detailed balan
e Eq.(8), we have
s,−1)
, whi
h when substituted in this equation gives
Π(s′) P(s′,−1→s′,1)+ P(s′,−1→s,−1) =Π(s′).
" #
s
X
s s = s′
Note that we have
ompleted the sum over with the term be
ause its
probability vanishes. The term in the bra
ket exhausts all possibilities for a move