Accurate nonadiabatic quantum dynamics on the cheap: making the most of mean field theory with master equations

0.72 MB·English

by  Aaron Kelly

Accurate nonadiabatic quantum dynamics on the cheap:
making the most of mean field theory with master equations
Aaron Kelly,1 Nora Brackbill,2 and Thomas E. Markland1,∗
1Department of Chemistry, Stanford University, Stanford, CA 94305, USA
2Department of Physics, Stanford University, Stanford, CA 94305, USA
(Dated: February 18, 2015)
In this article we show how Ehrenfest mean field theory can be made both a more accurate
andefficientmethodtotreatnonadiabaticquantumdynamicsbycombiningitwiththegeneralized
quantummasterequationframework. Theresultingmeanfieldgeneralizedquantummasterequation
5 (MF-GQME) approach is a non-perturbative and non-Markovian theory to treat open quantum
1 systems without any restrictions on the form of the Hamiltonian that it can be applied to. By
0 studying relaxation dynamics in a wide range of dynamical regimes, typical of charge and energy
2 transfer, we show that MF-GQME provides a much higher accuracy than a direct application of
b meanfieldtheory. Inaddition,theseincreasesinaccuracyareaccompaniedbycomputationalspeed-
e ups of between one and two orders of magnitude that become larger as the system becomes more
F nonadiabatic. This combination of quantum-classical theory and master equation techniques thus
makes it possible to obtain the accuracy of much more computationally expensive approaches at a
7
cost lower than even mean field dynamics, providing the ability to treat the quantum dynamics of
1
atomistic condensed phase systems for long times.

]
h PACSnumbers:
p
-
m I. INTRODUCTION ing further dynamical correlations between the subsys-
tem and the bath adds further accuracy, at the expense
e
h of assigning weights and phase factors to the trajectories
The exact treatment of real time nonadiabatic quan-
c [10–12]. This in turn adds many orders of magnitude to
tum dynamics in condensed phase chemical systems re-
s. mains a significant challenge that spurs the ongoing de- thenumberoftrajectories,whichgrowsrapidlywithsys-
c tem dimensionality and time, that must be generated in
i velopmentofapproximatemethodsthatareaccurate,ef-
s order to obtain converged properties.
ficient,andcantreatsystemswithawiderangeofdiffer-
y
h entformsofinteractions. Inparticularquantum-classical Sincemovingupthishierarchyrequiresordersofmag-
p (semiclassical)trajectorybasedmethodsofferahierarchy nitude more computational effort, only the lowest tiers
[ ofapproaches,derivedfromtheexactrealtimepathinte- are likely to be practical, both now and in the foresee-
gral formulation of quantum mechanics, that offer differ- able future, for nonadiabatic problems containing large
2
ent balances between accuracy and computational cost. quantumsubsystemsorwhereon-the-flytreatmentofthe
v
electronic states is required. At present, this means that
7 At the lowest tier of this hierarchy is Ehrenfest mean
4 field theory (MFT), which neglects all dynamical corre- one is typically limited to using MFT, or Tully’s fewest
1 switches surface hopping (FSSH) algorithm and its vari-
lations between between the quantum (subsystem) and
3 ants [13–17], or linearized path integral approaches [18–
classical(bath)degreesoffreedom[1]. Abovethislielin-
0 20].
earizedpathintegralapproachessuchasLSC-IVR[2,3],
.
1 FK-LPI [4], PBME [5], and LAND-map [6] which, al- Thegeneralizedquantummasterequation(GQME)of-
0 though still mean field in nature, capture some correla- fers an alternative way of exactly describing the nona-
5
tion. These methods offer higher accuracy than mean diabatic evolution of a quantum subsystem by formally
1
: field theory, at the expense of at least an order of mag- recasting the effects of the environment into a memory
v nitude more computational effort arising from the need kernel [21, 22]. In the condensed phase the environment
i
X to average over the mapping variables associated with comprisesofalargenumberofmodeswithabroadrange
the quantum subsystem degrees of freedom. Recently it of frequencies that couple to the quantum subsystem,
r
a has been shown that partially linearizing the propagator leading to a memory kernel that typically decays much
for the electronic degrees of freedom, giving rise to the more rapidly than the population relaxation time of the
partially linearized density matrix (PLDM) [7] approach subsystem. This separation in time-scales becomes more
andtheforward-backwardtrajectorysolution(FBTS)to pronounced as the system-bath coupling strength or the
the quantum-classical Liouvilleequation [8], allows more nonadiabaticity is increased. In such regimes, trajectory
dynamical correlation to be included at relatively little based approaches suffer from the rapid accumulation er-
additional cost to fully linearized methods. Introduc- rors in lower tier methods and a rapid rise in the com-
putational cost with time in higher tiered approaches.
Hence, using trajectory based approaches to calculate
the memory kernel, which is short lived compared to the
∗Electronicaddress: [email protected] subsystemrelaxationtime,andthenusingtheGQMEto
2
generate the subsystem dynamics, offers massive advan- The quantum-classical Liouville equation [26],
tages in terms of both accuracy and computational cost
∂
compared to a direct application of the trajectory based
ρˆ (X,t)=−iLρˆ (X,t), (3)
approach. Such a realization has previously been shown ∂t W W
to be highly effective in the case of higher-tier trajectory
describes the time evolution of the density matrix
based methods [23–25]. However, it also allows one to
ρˆ (X,t), which is a quantum mechanical operator that
W
make the most of the lower tier approaches by allowing
depends on the classical phase space variables. The
a more accurate treatment of strongly nonadiabatic and
quantum-classical Liouville (QCL) operator is
coupledproblemswhileretainingtheirexistingstrengths
in treating weakly coupled and adiabatic regimes. i 1
iL·= [Hˆ ,·]− ({Hˆ ,·}−{·,Hˆ }), (4)
Here we show that MFT can be combined with the (cid:126) W 2 W W
GQME formalism to yield a method that is as accu-
where [·,·] is the commutator, and {·,·} is the Poisson
rate as higher tiered trajectory-based techniques such as
bracketinthephasespaceoftheenvironmentalvariables.
FBTS and PLDM, but requires a much lower computa-
The subscript W refers to the partial Wigner transform
tional effort than MFT alone, without being limited to
overtheenvironmentaldegreesoffreedominthesystem.
any particular form of the Hamiltonian. We show how
The partial Wigner transform of the density operator, ρˆ,
to calculate the memory kernel of the GQME using dy-
is
namical trajectories obtained by solving the mean field
equationsofmotionforthesystem,andhowtousethese 1 (cid:90) Z Z
mean field kernels to subsequently generate the reduced ρˆW(R,P)= (2π(cid:126))Nb dZeiP·Z(cid:104)R− 2|ρˆ|R+ 2(cid:105). (5)
dynamicsofthequantumsubsystem. Finallywedemon-
strate that our method is more computationally efficient In order to arrive at MFT, one makes the approxima-
and more physically accurate for treating charge and en- tion that the density of the system can be written as a
ergy transfer regimes of the spin-boson model than a di- product of the subsystem and bath reduced densities at
rect application of MFT. alltimes. Thisassertsthattherearenodynamicalcorre-
lationsbetweenthesubsystemandthebath, asthetotal
density remains in a product state.
II. THEORY
ρˆ (X,t)=ρˆ (t)ρ (X,t), (6)
W s b
A. Ehrenfest Mean Field Theory
wherethebathRDMisρ (X,t)=Tr (ρˆ (X,t)). Using
b s W
the approximations in Eqs. (3) and (6), one obtains the
A particularly simple and instructive route to derive
Ehrenfest mean-field equations of motion [27]:
the Ehrenfest mean field equations of motion is to begin
with the quantum-classical Liouville equation [26] and ∂ρˆ (t) i (cid:20) (cid:18)(cid:90) (cid:19) (cid:21)
s = − Hˆ + dXHˆ (R(t))ρ (X,t) ,ρˆ (t) ,
neglect correlations in the system-bath dynamics. One ∂t (cid:126) s sb b s
begins by considering a system in which only a small
dR(t) P(t)
subset of the degrees of freedom behave quantum me- = , (7)
dt M
chanically, and are of interest. This set of degrees of
freedom is denoted as the quantum subsystem (or sim- dP = −∂(Hb,W(X(t))+Trs{Hˆsb(R(t))ρˆs(t)}).
ply as the subsystem), and the remainder of the system dt ∂R(t)
is referred to as the bath, and is assumed to behave es-
The mean field approximation for a subsystem observ-
sentially classically.
able, Oˆ (t), can be written as
The total Hamiltonian for the entire system is written s
as a sum of subsystem, bath, and coupling terms, (cid:90) (cid:16) (cid:17)
(cid:104)Oˆ (t)(cid:105)= dX Tr {Oˆ ρˆ (t)} ρ (X,t). (8)
Hˆ =Hˆ +Hˆ +Hˆ , (1) s s s s b
s b sb
where the subscripts s,b, and sb refer to the subsystem, MFT is a highly efficient approximation to quantum dy-
the bath, and the system-bath coupling, respectively. namics which, while obtaining accurate results in some
The time-evolution of the reduced density matrix of the regimes, fails when quantum effects in the bath are im-
subsystem is defined as portantandinregimeswithanonzerosubsystemenergy
(cid:90) bias.
ρˆ (t)=Tr (ρˆ(t))= dXρˆ (X,t), (2)
s b W
where ρˆ is the density operator for the entire system B. The Generalized Quantum Master Equation
and ρˆ is its Wigner transform, Tr indicates the par-
W b
tial trace taken over the bath degrees of freedom, X = An alternative approach to formulating the expecta-
(R,P) = (R ,R ,...,R ,P ,P ,...,P ), and N is the tionvalueofasubsystemobservableintermsoftheden-
1 2 Nb 1 2 Nb b
number of bath degrees of freedom. sity matrix for the full system, is to treat the reduced
3
dynamics of the subsystem via projection operator tech- initial state (that can also be calculated by mean field
niques. One begins with the exact quantum evolution of theory and other approximate methods). The free sub-
thefullsystem,governedbytheLiouville-vonNeumann system evolution prescribed by L is generally simple to
s
equation, simulate, and hence calculating the evolution of the sub-
system RDM reduces to the calculation of the memory
∂ i
ρˆ(t)=− [Hˆ,ρ(t)], (9) kernel, which encodes the deviation from free subsystem
∂t (cid:126) evolution. From Eq. (15), it is clear that changes in the
subsystem populations in the GQME are driven by the
and focuses only on the evolution of the subsystem de-
subsystem Liouville operator as well as the memory ker-
greesoffreedombyprojectingoutthedegreesoffreedom
nel and thus the lifetime of the population dynamics is
ofthebath. Thiscanbeaccomplishedbyapplyingapro-
typically longer than that of the memory kernel. In con-
jection operator, P,
densed phase systems, where the environment spans a
(cid:90)
wide spectral bandwidth, the memory kernel is expected
P =ρˆeq⊗Tr (·)=ρeq dX(·), (10)
b b b,W to decay much more quickly than the population relax-
ation time.
and its compliment, Q=1−P. The general form for the memory kernel, given above,
TheformofthecouplingpartoftheHamiltonian, Hˆsb is not straightforward to evaluate since it explicitly de-
is chosen to be pends on the projection operator. However, there are a
variety of ways that it can be constructed from simula-
Hˆ =Sˆ⊗Λˆ, (11)
sb tionsoftheunprojecteddynamicsofthefullsystem. We
where Sˆ is a pure subsystem operator, and Λˆ is a pure will limit our discussion of this procedure to a particular
method introduced by Shi and Geva, who showed that
bath operator. We use this factorization for simplicity
the full memory kernel, K(τ), can be written in terms of
and, since any coupling Hamiltonian can be written as a
a set of partial memory kernels [23, 28, 29],
sum of such operators, there is no loss of generality [28].
WealsowritethebathpartofthecouplingHamiltonian, (cid:90) τ
Λˆ, such that its equilibrium thermal average vanishes, K(τ) = K (τ)+i dτ(cid:48)K (τ −τ(cid:48))K (τ), (17)
1 1 2
0
(cid:104)Λˆ(cid:105) =Tr [Λˆρˆeq]=0, (12) (cid:90) τ
eq b b K (τ) = K (τ)+i dτ(cid:48)K (τ −τ(cid:48))K (τ), (18)
2 3 3 2
In the problems that we consider in Sec. III such a con- 0
dition is naturally satisfied, but in general this condition where the partial memory kernels are given by
can be enforced by redefining Λˆ relative to its thermal
average [29]. K (τ) = Tr {L e−iLτL ρˆeq}, (19)
1 b sb sb b
We restrict our considerations to initial states of the K (τ) = Tr {e−iLτL ρˆeq}. (20)
Feynman-Vernon type, where the initial density of full 3 b sb b
system can be factorized in the following manner,
By combining Eq. (15), and Eqs. (17 - 20) the subsys-
ρˆ(t=0)=ρˆ (0)⊗ρˆeq, (13) tem RDM can thus be exactly evolved using projection-
s b free input. Exact numerical evaluations of these expres-
where sionshaverecentlybeencarriedoutusingtechniquessuch
as QUAPI, real-time quantum Monte Carlo, and ML-
exp(−βHˆ ) MCTDH [29–32].
ρˆeq = b (14)
b Tr [exp(−βHˆ )]
b b
is the density operator for the isolated bath in thermal C. Mean Field Evaluation of the Memory Kernel
equilibrium. Under these conditions the exact time evo-
lution of the subsystem RDM is given by the Nakajima- In practice, evolving the subsystem RDM using the
Zwanzig GQME [21, 22] GQMEformalismforanarbitrarysystemisnolesscum-
bersome than solving Eq. (9). One way to proceed, that
d (cid:90) t
ρˆ (t)=−iL ρˆ (t)− dτK(τ)ρˆ (t−τ), (15) is valid in the weak coupling limit, is to factorize the
dt s s s 0 s memorykernelintosubsystemandbathpartswhichcan
be evaluated separately, leading to the well established
where the subsystem Liouville operator, L , is given by
L = 1[Hˆ ,·], and the memory kernel, K, issgiven by Bloch-Redfield theory [29, 33]. In this limit K2 and K3
s (cid:126) s vanish, and the memory kernel can be evaluated using
K(τ)=Tr {L exp(−iQLτ)QL ρˆeq} (16) equilibrium correlation functions of the isolated bath. In
b sb sb b
contrast, here we make no such simplifying assumption
Theassumptionofaninitiallyuncorrelatedstateisnot and instead evaluate the kernels using MFT.
necessary, and relaxing it would introduce an inhomoge- Since K and K do not contain any projected input,
1 3
neous term to the GQME which depends on the chosen they can be simulated directly and then used to obtain
4
K by solving Eqs. (17) and (18). The matrix elements 2. K is generated from K by an iterative solution
2 3
of K and K are obtained by projecting each quantity to Eq. (18), using K itself as an initial guess for
1 3 3
onto a basis which spans the subsystem Hilbert space, K . This iterative procedure typically converges
2
very quickly, and often requires only a few tens of
(cid:68) (cid:69)
(K1)αα(cid:48)ββ(cid:48)(τ) = Sαµ(cid:48)(τ)Λˆβµ(cid:48)(cid:48)µα(cid:48)(τ)Sµβ(0)Λˆ(0) iterations.
eq
(cid:68) (cid:69)
− S (τ)Λˆβ(cid:48)µ(cid:48)(τ)S (0)Λˆ(0) 3. K1andK2areusedasinputtoobtainthefullmem-
µ(cid:48)α(cid:48) αµ µβ eq ory kernel K by numerical integration of Eqs. (17)
+(cid:68)Λˆ(0)S (0)Λˆµµ(cid:48)(τ)S (τ)(cid:69) and (18).
β(cid:48)µ αβ µ(cid:48)α(cid:48)
eq
(cid:68) (cid:69) 4. Using the full memory kernel, the evolution of the
− Λˆ(0)Sβ(cid:48)µ(0)Λˆµµα(cid:48)β(cid:48)(τ)Sαµ(cid:48)(τ) , subsystemdensityisgeneratedbydirectnumerical
eq
integration of the GQME using Eq. (15).
(21)
In our calculations the kernel elements were calculated
(cid:68) (cid:69) for the specified time τ, and set to zero for all t>τ; no
(K ) (τ) = (ˆ1 )β(cid:48)α(cid:48)(τ)S (0)Λˆ(0)
3 αα(cid:48)ββ(cid:48) b αµ µβ smoothing of the kernel data was performed for t < τ.
eq
(cid:68) (cid:69) Using the MF-GQME approach one can propagate the
− S (0)Λˆ(0)(ˆ1 )µα(cid:48)(τ) ,(22)
β(cid:48)µ b αβ subsystem RDM for arbitrarily long times using only
eq
short-time information obtained from mean field trajec-
where α,α(cid:48),β, and β(cid:48) refer to subsystem states, and ˆ1 tories.
b
is the unit operator for the bath. In the above two ex-
pressions the Einstein summation convention is used.
The matrix elements of the partial memory kernels, III. RESULTS AND DISCUSSION
K andK , containcorrelationfunctionsofthefollowing
1 3
form, In order to assess the accuracy and efficiency of our
(cid:16) MF-GQME approach, we performed simulations of the
(cid:104)Sˆ(0)Λˆ(0)Γˆβα(cid:48)βα(cid:48)(τ)(cid:105)eq = Tr ρˆebqHˆsb|β(cid:105)(cid:104)β(cid:48)| spin-boson model. Despite its apparent simplicity, this
(cid:17) system is a prototypical model for the study of quan-
×eiLτ/(cid:126)Γˆ|α(cid:48)(cid:105)(cid:104)α| , (23)
tumtransportandrelaxationprocessesinthecondensed
phase [34, 35], and remains a challenging test to ap-
whereΓˆ isabathoperator(ˆ1b orΛˆ). Definingoperators, proximate methods. Since it is now possible to gener-
A = Hˆ (ˆ1 ⊗|β(cid:105)(cid:104)β(cid:48)|) and B = ˆ1 ⊗|α(cid:48)(cid:105)(cid:104)α| , or B = ate numerically exact results in many of the parameter
sb b b
Hˆ (ˆ1 ⊗|α(cid:48)(cid:105)(cid:104)α|), expression (23) takes on the general regimes of the spin-boson model, it provides an ideal
sb b
form for a quantum time correlation function, benchmark test case for the accuracy and efficiency of
approximate nonadiabatic dynamics approaches. In par-
(cid:104)Sˆ(0)Λˆ(0)Γˆβ(cid:48)α(cid:48)(τ)(cid:105) =Tr(ρˆeqAB(τ)). (24) ticular we compare our MF-GQME approach to a direct
αβ eq b
MFT treatment, as well as to the recently introduced
where the equilibrium density corresponds to that of the
FBTSmethod[8,9]whichhasbeenshowntooutperform
isolated bath.
fullylinearizedmethodsatmarginalextracomputational
Working in the coordinate representation of the bath
cost.
degreesoffreedom,andmakinguseofthepartialWigner
Thespin-bosonHamiltoniancanbewritteninthesub-
transform, Eq. (24) can be rewritten as
system basis as
(cid:90)
Tr(ρˆebqAB(τ))=Trs dX[ρˆebqA]W (X,0)BW(X,τ)(.25) Hˆ =(cid:15)σˆ +∆σˆ + Pˆ2 +(cid:88)(cid:18)1M ω2Rˆ2−c Rˆ σˆ (cid:19),
z x 2M 2 j j j j j z
The calculation of the memory kernel then amounts j
(26)
to the evaluation of the above expression, which can be
where σˆ and σˆ are Pauli spin matrices. This Hamil-
performed by a hybrid Monte Carlo / molecular dynam- x z
tonian describes a two level quantum system with ener-
ics algorithm where (i) initial conditions are sampled
from [ρˆeqA] (X,0), and (ii) the system is propagated getic bias 2(cid:15), and electronic coupling (tunneling) matrix
b W element ∆, that is bi-linearly coupled to a bath of inde-
intimefromtheseinitialconditionsusingMFTtoevalu-
pendent harmonic oscillators. In the spin-boson model
ate B (X,τ). The full memory kernel can then be con-
W both subsystem states are coupled to the same bath in
structed,andthesubsystemRDMpropagatedasfollows:
ananti-correlatedfashion. Incontrasttoproblemswhere
1. Mean-fieldtrajectoriesareusedtoobtainthecorre- the subsystem states are coupled to independent uncor-
lationfunctionsnecessarytoformK andK using related baths, such as the Frenkel exciton model, this
1 3
Eqs. (21) and (22). See Appendix A for more de- form of the coupling presents a more difficult challenge
tails on the calculation of these quantities for the for mean-field type methods to describe [9, 36–41]. The
model studied here. greater challenge of treating anti-correlated baths is due
5
1 0.5 Re K Im K
1211 0.3 1211
ξ = 0.1
0
0
0.5
-0.5
-0.3
0 Re K Im K
1212 1212
0 0.5
> -0.5 -2
0

(t) 0 5 10 15 20 -4
0.3 Re K 0.3 Im K
σz 1 1221 1221
< ξ = 0.4 0 0
0.5
-0.3 -0.3

0 0.5 Re K1222 0.5 Im K1222

0 0

-0.5

-0.5 -0.5

-1 0 2 4 0 2 4
0 2 4 6 8 10 t ∆ t ∆

t ∆

Figure2: MatrixelementsofthememorykerneloftheGQME

forforωc =2∆,(cid:15)=∆,β =5∆−1,ξ=0.4,δ=0.02∆−1,and
N =2x104. In each panel the MFT results are shown as
traj
Figure 1: Evolution of the subsystem population difference blue lines, and the exact QUAPI results are the black lines.
in the biased, nonadiabatic regime with increasing system- Notethedifferenty-axisscalesineachpanel,duetothevary-
bath coupling strength; ωc = 2∆, (cid:15) = ∆, β = 5∆−1. In ing magnitudes of each element.
eachpaneltheexactresultsareshowninsolidblackdots,the
dotted blue lines are direct MFT results, the solid red lines
are FBTS results, and the solid blue lines are the results of
our MF-GQME approach. integral (QUAPI) algorithm [42–44]. The performance
of trajectory based approaches degrades as the system
becomes more nonadiabatic and the system-bath cou-
to the greater difference between the mean force and pling is increased, a failure that is particularly evident
those on the individual diabatic surfaces compared to for systems with nonzero energetic bias. In this regime
independent, uncorrelated, baths. direct MFT is over-coherent and only captures the exact
The interaction between the system and the bath can QUAPIresultsatveryshorttimes(τ∆<1.5). Thelong
be fully characterized by the spectral density, J(ω), time population difference is close to zero, which is a no-
which determines the strength of the interactions be- toriousfeatureofMFTthatarisesduetothedeviationof
tween the subsystem and bath, which we chose to be themeanforcefromtheforcesontheindividualsurfaces,
of Ohmic form, resulting in an accumulation of error in the subsystem
populationdistributionastimeprogresses. FBTS,which
π
J(ω) = ξωe−ω/ωc. (27) is also mean field in nature but includes more dynamical
2
correlation between the subsystem and bath, much more
The Kondo parameter, ξ, controls the strength of the accurately captures the long time populations at a mod-
coupling between subsystem and the bath, and the cut- est increase in the number of trajectories of by a factor
off frequency ω sets the primary time-scale for the bath of ≈ 5 in this regime. In contrast to MFT, the FBTS
c
evolution. The quantum subsystem was initialized in di- underestimates the oscillatory nature of the population
abatic state 1, and the bath was sampled from its (iso- decay.
lated) equilibrium distribution. In our calculations 400 By using MFT to approximate the memory kernel,
bath modes were used to represent the continuous spec- our MF-GQME approach produces populations dynam-
tral density which, for all regimes and approaches em- ics that are in perfect quantitative agreement with the
ployed, gave results converged to graphical accuracy. exact quantum results. Figure 2 shows the MFT mem-
Figure 1 compares the subsystem evolution as the ory kernel elements that give rise to the dynamics in the
system-bath coupling is increased, for a system with bottompanelofFig. 1,comparedwiththeexactQUAPI
an energetic bias ((cid:15) (cid:54)= 0) in the nonadiabatic regime results. There are four nonzero, linearly independent el-
(ωc > 1). Numerically exact results were generated us- ements of K for a two-level system (due to the form of
∆
ing our own implementation of the quasi-adiabatic path- the spin-boson Hamiltonian, K = 0). As expected,
ααββ(cid:48)
6
1 1
ω = ∆
c
0.75
0.5
0.5
0.25
0
> 0
1 )
> ωc = 2.5∆ (tz -0.251 ωc = ∆
σ τ ∆ = 0.1
(t) < ττ ∆∆ == 00..35
σ z 0 0.5 ττ ∆∆ == 00..79

< 0 τM ∆F T= 1.5

1

ω = 7.5∆
c -0.5 ω = 7.5 ∆
c

0 0 5 10 15
t ∆

0 5 10 15 Figure 4: Evolution of the subsystem population difference
t ∆ generated from the MF-GQME approach using memory ker-
nels of varying length; (cid:15)=∆, β =5∆−1, ξ=0.1.
Figure3: Evolutionofthesubsystempopulationdifferencein
the intermediate coupling regime with increasing bath nona-
diabaticity; (cid:15) = ∆, β = 5∆−1, ξ = 0.1. In each panel the performsbetterthandirectMFTinreproducingthelong
exact results are the solid black dots, the dotted blue lines timelimit,albeitatanincreasedcostofanorderofmag-
are direct MFT results, the solid red lines are FBTS results,
nitude more trajectories, but again gives results which
andthesolidbluelinesaretheresultsofourMF-GQMEap-
are increasingly overdamped as ω is increased. This
proach. c
overdamping is consistent with that observed in previ-
ous FBTS and PLDM results in similar regimes [9, 40].
MFT fails to correctly capture the long time behavior; The MF-GQME results are in quantitative agreement
exhibiting spurious oscillations when τ∆ > 2. This is at low nonadiabaticity, and even at the highest nonadi-
consistent with the fact that the direct MFT treatment abaticity exhibit only a very subtle phase shift relative
ofthepopulationdynamicsinFig. 1alsobeginstoshow to the exact results. Again this reflects that the mem-
marked errors at those times. Despite the deviations of ory kernel decays rapidly in these nonadiabatic regimes,
MFT from the exact kernels, the population dynamics and hence retaining a memory kernel of length τ∆=1.5
generated using a memory kernel of length t∆ = 1.5 is sufficient to generate the results shown in Fig. 3.
(shown in the bottom panel of Fig. 1) are in excellent Figure 4 shows the convergence of the population dy-
agreement with the exact results. Using memory times namics obtained from the MF-GQME approach when
longerthan τ∆=1.5includes thespuriouslong timeos- different amounts of time, τ∆, are included in the mem-
cillations but only introduces very minor changes to the ory kernel for the regimes shown in the top and bottom
population dynamics. As MF-GQME only requires the panels of Fig. 3. In both cases the convergence is essen-
generation of very short trajectories to obtain the mem- tially monotonic as the length of memory in the kernel
ory kernel, the entire population decay, which occurs in is increased. The MF-GQME results are generally bet-
a time of approximately 15∆−1, can be obtained at a terthandirectMFTevenforveryshortmemorykernels.
cost 10 times cheaper than a standard mean-field calcu- When an insufficient length of time is used in the ker-
lation of the same observable and ≈ 50 times cheaper nel,theerroraccumulatedinthepropagationofthesub-
than FBTS. system RDM in MF-GQME manifests in the observed
The effect of increasing the characteristic frequency of population dynamics more prominently at longer times.
the bath, ω , which pushes the system into an increas- For the more adiabatic regime shown (top panel) the
c
ingly nonadiabatic regime, is displayed in Fig. 3. In the dissipation induced by the bath is well captured when
adiabaticlimitMFTisaccurate,butasthenonadiabatic- τ∆=0.7,whereaswhenthesystembecomesmorenona-
ity is increased the direct MFT results begin to deviate diabatic (bottom panel) this is decreased further with
from the exact solution at progressively shorter times, convergence obtained around τ∆ = 0.5, which is 30
and by ω = 7.5∆ it is unable to even capture the first times shorter than the population decay time. This
c
minimum in the coherent decay correctly. FBTS again again highlights the complementarity of using MFT and
7
1
ω = 5∆, ξ = 0.2 1
c
0.5
0.5
q
e
>
0 z 0
σ

-0.5 < -0.5
1
-1
ω = 5∆, ξ = 0.5
> c
-2

t) 0.5
(
σz ∆ )
< / -3
0 k

(
1 n
l
ω = 10∆, ξ = 0.4 -4
c
0.5
-2 0 2
2 ε / λ
0
Figure 6: Marcus electron transfer regime, with ω = 10∆,
c
ξ = 1.0, β∆ = 0.05, and λ = 2ξω is the reorganization
c
0 2 4 6 8 10 energy. Thetoppanelshowstheequilibriumsubsystempop-
t ∆ ulationdifferenceafterrelaxationhasoccurredwhilethebot-
tom panel shows the electron transfer rate as a function of
drivingforceobtainedfromexponentialfitstothepopulation
Figure 5: Evolution of the subsystem population difference
decay. In the top panel the solid black line is the Boltzmann
at zero temperature in the unbiased case, (cid:15) = 0. In each
distribution, in the bottom panel the solid black line is the
panel the exact ML-MCTDH results from Ref [45] are the
Marcusrate,thedottedbluelinesareMFTresults,thesolid
solidblackdots,thedottedbluelinesaredirectMFTresults,
red circles are FBTS results, and the solid blue squares are
the solid red lines are FBTS results, and the solid blue lines
the results of our MF-GQME approach.
are the results of our MF-GQME approach.
FBTS performs worse than MFT in these regimes and
other trajectory based approaches in conjunction with again exhibits over-damped behavior resulting in slower
the GQME framework; moving further into the nonadia- relaxation. TheMF-GQMEresultsareinexcellentagree-
baticregime,whichleadstoafasterbreakdownofMFT, ment, with some very mild underdamping present in the
requireslesstotaltimetobeincludedintheGQMEmem- case where the system-bath coupling is strong (middle
ory kernel. panel).
In contrast to the spin-boson model at finite tempera- In addition to the significant increase in accuracy af-
ture,thezero-temperaturelimitisconsideredtobemore fordedbyMF-GQMEoverFBTSanddirectMFT,thisis
challenging, as nuclear quantum effects in the bath be- accompanied by a significant increase in efficiency. As in
come more prominent; a classical treatment of the bath Fig. 4,thememorykerneldecaytimesinFig. 5required
distributionwouldpredictonlyasingleallowedbathini- to obtain the converged population dynamics are sub-
tial configuration. In all the results shown here, the stantially shorter than the population decay times they
Wigner sampling of the bath ensures that zero-point en- generate, with memory kernel of lengths τ∆ = 0.2, 0.2
ergyisincludedexactlyintheinitialcondition,although and 0.1 required to obtain the results in the top, mid-
this is not guaranteed to be preserved by the approx- dle and bottom panels respectively. This again reflects
imate evolution prescribed by MFT or FBTS. In Fig. that the memory kernel becomes shorter as the system
5 we compare to exact multi-layer multi-configurational becomesmorenon-adiabatic. Duetothemuchfasterde-
time-dependent Hartree (ML-MCTDH) results at zero- cayofthememorykernelthanthepopulationsdynamics
temperature in the nonadiabatic regime [45]. As seen in (which occurs in approximately 5∆−1), the MF-GQME
Fig. 1 and Fig. 3 increasing the nonadibaticity, ω /∆, results in Fig. 4 are between 25 and 50 times cheaper to
c
or system-bath coupling, ξ, degrades the performance of generatethandirectMFTdynamicsandupto500times
MFT and FBTS at progressively shorter times although cheaper than FBTS.
due to the lack of energetic bias the long time limits are One of the most well-known difficulties of mean field
in, probably fortuitously, good agreement. Suprisingly, theory occurs in the Marcus electron transfer regime
8
[46]. Whilemostapproximatedynamicsapproaches, like and only the short-lived memory kernel term is approx-
FSSH and MFT, are capable of qualitatively capturing imated by MFT, which is accurate at short-times. The
the famous rate turnover as a function of driving force complementary nature of the combination of MFT with
[16, 25, 47, 48], as shown in the bottom panel of Fig. 6, the GQME is amplified as the dynamics become more
MFTfailstocorrectlydescribethedonor-acceptorprod- strongly nonadiabatic, as the memory kernel becomes
uctratiosfortheprocess. MuchlikethePBMEapproach increasingly short-lived compared to the subsystem pop-
[25], the FBTS method qualitatively captures the rate ulation relaxation time. While MFT can fail at long-
turnoverbuthasaslightlyasymmetricshape. Also, per- times when used directly, the timescale separation be-
haps surprisingly given their mean field treatment of the tween the population dynamics and the memory ker-
system-bath interactions, both PBME and FBTS give nel decay in these cases ensures that the dynamics can
the correct equilibrium distribution at long times. The still be captured accurately using the MF-GQME. The
ability of FBTS to qualitatively and semiquantitatively MF-GQMEapproachisthereforeexpectedtobeefficient
captureMarcusturnoverisconsistentwithrecentPLDM enoughbeusedtostudynonequilibriumrelaxationprob-
results for a similar model of Marcus electron transfer lems in complex condensed phase systems at a very low
that was simulated using a flux-side correlation function computational cost, and it can be applied to any form of
approach[49]. TheagreementbetweenFBTSandPLDM the system, bath or coupling between them, i.e. it is in
is expected as the methods are identical when the sub- no way limited to the linear coupling or harmonic bath
system Hamiltonian is chosen to be traceless [9]. They invoked in the spin-boson model studied here.
thus offer a similar tradeoff between cost and accuracy Although MF-GQME offers probably the highest pos-
in the semiclassical hierarchy. sible accuracy for the lowest possible cost, other simi-
While FBTS performs well by including some dynam- lar calculations can be carried out using other dynam-
ical correlations between the system and bath, the com- ics methods such as FSSH, linearized and partially lin-
pleteneglectofthesecorrelationsindirectMFTleadsto earized approaches, and higher tier methods [2, 5, 7,
poorperformance. Thisisstronglypronouncedinbiased 8, 17]. Indeed, at the opposite end of the hierarchy,
regimes(thosewithlargeelectronicdrivingforces)where very computationally demanding approaches such as the
the reaction rates are underestimated and the long time momentum-jumpsolutiontotheQCLE[11]canbemade
population difference is too small. The MF-GQME ap- tractable when combined within the GQME framework
proach resolves these issues, giving rise to quantitative [25], greatly expanding its regime of applicability. In
agreement with the Marcus prediction for the rates and addition, the method for generating the memory ker-
theBoltzmanndistributionofthelongtimepopulations. nel from system-dependent bath correlation functions
In order to obtain the MF-GQME results in this regime adopted here represents just one possible way obtaining
τ∆ = 0.5 time units of data were typically included in K from semiclassical trajectory-based simulation meth-
the memory kernels, and the population decay occurs on ods, and future work will explore these issues.
atimescaleofafewhundredtimeunits. Thisallowedfor
a ≈200 fold efficiency gain compared to direct MFT, in
addition to a substantially improved accuracy.
V. ACKNOWLEDGEMENTS
IV. CONCLUSIONS The authors would like to thank Tim Berkelbach,
Will Pfalzgraff and Hans Andersen for helpful comments
Here we have shown that utilizing MFT within the and a thorough reading of this manuscript. T.E.M. ac-
GQME framework gives rise to nonadiabatic relaxation knowledges support from a Terman fellowship, an Al-
dynamics that are highly accurate across a wide range fred P. Sloan Research fellowship, a Hellman Faculty
of physical regimes. This is accompanied by a computa- ScholarFundfellowshipandStanfordUniversitystart-up
tional savings of one or two orders of magnitude com- funds. A.K.acknowledgesapostdoctoralfellowshipfrom
pared with direct MFT calculations, and up to three the Stanford Center for Molecular Analysis and Design.
orders of magnitude over FBTS. This success can be N.B. was supported by the National Science Foundation
rationalized based on the fact that the subsystem Li- GraduateResearchFellowshipProgramunderGrantNo.
ouville operator is treated exactly within the GQME DGE-114747.
[1] A. D. McLachlan Mol. Phys. 8, 39 (1964). [5] H.Kim,A.Nassimi,andR.Kapral.,J.Chem.Phys.129,
[2] X. Sun, H. B. Wang, and W. H. Miller, J. Chem. Phys. 084102 (2008).
109, 7064 (1998). [6] S.BonellaandD.F.Coker,J.Chem.Phys.122,194102
[3] Q. Shi and E. Geva, J. Chem. Phys. 118, 8173 (2003). (2005).
[4] J. A. Poulsen, G. Nyman, and P. J. Rossky, J. Chem. [7] P. Huo and D. F. Coker, J. Chem. Phys. 135, 201101
Phys. 119, 12179 (2003). (2011).
9
[8] C.-Y.HsiehandR.Kapral.,J.Chem.Phys.137,22A507 [44] N.Makri,andD.E.Makarov,J.Chem.Phys.102,4611
(2012). (1994).
[9] C.-Y.HsiehandR.Kapral.,J.Chem.Phys.138,134110 [45] H.Wang,andM.Thoss,NewJ.Phys.10,115005(2008).
(2013). [46] R.A. Marcus Annu. Rev. Phys. Chem. 15, 155 (1964).
[10] D. MacKernan, G. Ciccotti, and R. Kapral., J. Phys.: [47] W. Xie, S. Bai, L. Zhu, and Q. Shi, J. Phys. Chem. A,
Condens. Matter 14, 9069 (2002). jp400462f (2013).
[11] D. MacKernan, R. Kapral.,and G. Ciccotti, J. Phys. [48] C. A. Schwerdtfeger, A. V. Soudackov,and S. Hammes-
Chem. B 112, 424 (2008). Schiffer, J. Chem. Phys. 140, 034113 (2012).
[12] E. Dunkel, S. Bonella, and D. F. Coker, J. Chem. Phys. [49] P.Huo,T.F.MillerIIIandD.F.Coker,J.Chem.Phys.
129, 114106 (2008). 139, 151103 (2012).
[13] J. C. Tully,and R. K. Preston, J. Chem. Phys. 55, 562 [50] K. Imre, E. O¨zizmir, M. Rosenbaum, and P. F. Zwiefel,
(1971). J. Math. Phys. 5, 1097 (1967).
[14] J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
[15] E.R.Bittner,andP.J.Rossky,J.Chem.Phys.103,8130
(1995).
VI. APPENDIX
[16] B. R. Landry, and J. E. Subotnik J. Chem. Phys. 137,
22A513 (2012).
[17] B. R. Landry, M. J. Falk, and J. E. Subotnik J. Chem. A. Explicit expressions for computing K1 and K3
Phys. 139, 211101 (2013). for the spin-boson model
[18] W. H. Miller, J. Phys. Chem. A 105, 2942 (2001).
[19] M.Thoss,andH.WangAnn.Rev.Phys.Chem.55,299 AccordingtotheprocedureoutlinedinSec. (IIC),the
(2004).
integrationinEq. (25)iscarriedoutoverthephasespace
[20] R. Kapral, Annu. Rev. Phys. Chem. 57, 1239 (2006).
of the bath by Monte Carlo sampling initial conditions
[21] S. Nakajima, Prog. Theor. Phys. 20, 948 (1958).
from [ρˆeqA] (X,0) and generating dynamical trajecto-
[22] R. Zwanzig, J. Chem. Phys. 33, 1338 (1960). b W
ries to evaluate B (X,τ).
[23] Q. Shi and E. Geva, J. Chem. Phys. 120, 10647 (2004). W
[24] Q. Shi and E. Geva, J. Chem. Phys. 121, 3393 (2004). If A is independent of momentum and linear in the
[25] A. Kelly and T. E. Markland, J. Chem. Phys. 139, coordinatesofthebath, theWignertransformoftheop-
014104 (2013). erator product is [50]
[26] R. Kapral and G. Ciccotti, J. Chem. Phys. 110, 8919
(1999). i ∂A ∂ρeq
[27] R. Grunwald, A. Kelly, and R. Kapral, in Energy [Aρˆebq]∗W =[ρˆebqA]W =AWρeb,qW + 2(cid:126) ∂RW · ∂bP,W,
Transfer Dynamics in Biomaterial Systems, edited by (28)
I. Burghardt (Springer, Berlin, 2009), pp. 383–413.
[28] M.-L.Zhang,B.J.KaandE.Geva,J.Chem.Phys.125,
which can be evaluated analytically for the harmonic
044106 (2006).
bath employed in this study.
[29] Q. Shi and E. Geva, J. Chem. Phys. 119, 12063 (2003). In the spin-boson model, Λˆ = −(cid:80) c Rˆ and it’s
[30] L. Mu¨hlbacher, and E. Rabani, Phys. Rev. Lett. 100, j j j
(cid:80)
176403 (2008). Wigner transform is ΛW = − jcjRj. The Wigner
[31] G. Cohen and E. Rabani, Phys. Rev. B 84, 075150 transformoftheequilibriumdensityfortheisolatedbath
(2011). is
[32] E. Y. Wilner, H. Wang, M. Thoss,and E. Rabani, Phys.
Rev. B 90, 115145 (2014). ρeq (X)=(cid:89)tanhβωj/2
[33] A. G. Redfield, IBM J. Res. Dev. 1, 19 (1957). b,W π
[34] A. Leggett, S. Chakravarty, A. Dorsey, M. Fisher, A. j
[35] UG.arWg,eaisnsd, QR.uaZnwteurmgerD,iRsseivp.atMivoedS.yPshteyms.s59(W, 1or(l1d98S7c)i.en- ×exp(cid:34)−2tanhβωj/2(cid:32) P2 + Mjωj2R2(cid:33)(cid:35). (29)
ω 2M 2 j
tific, Singapore, 1992). j j
[36] G. Tao and W. H. Miller, J. Phys. Chem. Lett. 1, 891
(2010). In practice, our initial conditions for the bath degrees
[37] A. Kelly and Y. M. Rhee, J. Phys. Chem. Lett. 2, 808 offreedomweresampledfromtheinitialbathdensity(by
(2011). taking out a factor of the bath density from Eq. (28))
[38] T. C. Berkelbach, D. R. Reichman,and T. E. Markland, and the trajectories were then re-weighted using the re-
J. Chem. Phys. 136, 034113 (2012).
maining term in Eq. (28).
[39] T. C. Berkelbach, T. E. Markland,and D. R. Reichman,
The subsystem operator in the system-bath coupling
J. Chem. Phys. 136, 084104 (2012). part of the Hamiltonian is Sˆ = σˆ , and the matrix ele-
[40] P. Huo and D. F. Coker, J. Chem. Phys. 137, 22A535 z
ments of Sˆ in the subsystem basis are given by
(2012).
[41] P. Huo and D. F. Coker, J. Chem. Phys. 133, 184108
(cid:104)m|Sˆ|n(cid:105)=(cid:104)m|σˆ |n(cid:105)=S δ =(−1)n+1δ . (30)
(2011). z n mn mn
[42] D. E. Makarov, and N. Makri, Chem. Phys. Lett. 221,
482 (1994). where,mandnareeigenstatesoftheisolatedsubsystem,
[43] N.Makri,andD.E.Makarov,J.Chem.Phys.102,4600 i.e. σˆz|1(cid:105)=|1(cid:105) and σˆz|2(cid:105)=−|2(cid:105).
(1994). The elements of K and K are then given by
1 3
10
(cid:90) (cid:34) i ∂Λ ∂ρeq (cid:35)
(K ) (τ) = ((−1)β+α+(−1)β+α(cid:48)+1) dX Λ ρeq − W · b,W (X,0)Λ (X,τ)ρββ(cid:48)(0)ρα(cid:48)α(τ)
1 αα(cid:48)ββ(cid:48) W b,W 2(cid:126) ∂R ∂P W s s
(cid:90) (cid:34) i ∂Λ ∂ρeq (cid:35)
+((−1)β(cid:48)+α+1+(−1)β(cid:48)+α(cid:48)) dX Λ ρeq + W · b,W (X,0)Λ (X,τ)ρββ(cid:48)(0)ρα(cid:48)α(τ),
W b,W 2(cid:126) ∂R ∂P W s s
(31)
and
(cid:90) (cid:34) i ∂Λ ∂ρeq (cid:35)
(K ) (τ) = (−1)β+1 dX Λ ρeq − W · b,W (X,0)ρββ(cid:48)(0)ρα(cid:48)α(τ)
3 αα(cid:48)ββ(cid:48) W b,W 2(cid:126) ∂R ∂P s s
(cid:90) (cid:34) i ∂Λ ∂ρeq (cid:35)
+(−1)β(cid:48) dX Λ ρeq + W · b,W (X,0)ρββ(cid:48)(0)ρα(cid:48)α(τ), (32)
W b,W 2(cid:126) ∂R ∂P s s
where ραα(cid:48) =c c∗ is the subsystem RDM, and c and c are the coefficients of the subsystem wavefunction.
s α α(cid:48) α α(cid:48)