Non-Markovian quantum dynamics: Correlated projection superoperators and Hilbert space averaging

0.29 MB·English

by  Heinz-Peter Breuer

Non-Markovian quantum dynamics: Correlated projection superoperators and
Hilbert space averaging
Heinz-Peter Breuer,1,∗ Jochen Gemmer,2,† and Mathias Michel3,‡
1Physikalisches Institut, Universit¨at Freiburg, Hermann-Herder-Str. 3, D-79104 Freiburg, Germany
2Fachbereich Physik, Universit¨at Osnabru¨ck, Barbarastr. 7, D-49069 Osnabru¨ck, Germany
3Institut fu¨r Theoretische Physik, Universit¨at Stuttgart,
Pfaffenwaldring 57, D-70550 Stuttgart, Germany
(Dated: February 9, 2008)
Thetime-convolutionless(TCL)projectionoperatortechniqueallowsasystematicanalysisofthe
non-Markovianquantumdynamicsofopensystems. Weintroduceaclassofprojection superopera-
torswhichproject thestatesofthetotalsystemontocertaincorrelated system-environmentstates.
It is shown that the application of the TCL technique to this class of correlated superoperators
6
enables the non-perturbative treatment of the dynamics of system-environment models for which
0
the standard approach fails in any finite order of the coupling strength. We demonstrate further
0
that the correlated superoperators correspond to the idea of a best guess of conditional quantum
2
expectations which is determined by a suitable Hilbert space average. The general approach is

n illustratedbymeansofthemodelofaspinwhichinteractsthroughrandomlydistributedcouplings
a with a finite reservoir consisting of two energy bands. Extensive numerical simulations of the full
J
Schr¨odingerequation of the model reveal thepower and efficiency of themethod.

7
2 PACSnumbers: 03.65.Yz,05.70.Ln,05.30.-d

2
v I. INTRODUCTION pletely positive quantum dynamical semigroups and the
7 corresponding Markovian master equations in Lindblad
7
form [11, 12, 13, 14]. However,if the physical conditions
1 Realistic quantum mechanical systems are influenced
underlying the Markov approximation are violated one
0 through the coupling to an environment which contains
1 a large number of mostly uncontrollable degrees of free- has to cope with strong non-perturbative and memory
5 effects and the theoretical and mathematical treatment
dom. The unavoidable interaction of an open quantum
0 of the reduced system dynamics is typically much more
systems with its environment gives rise to the mecha-
/ involved.
h nisms of damping and dissipation, and to a strong and
p often rapid loss of quantum coherence. Applications of A systematic approach to non-Markovian dynamics is
t- the theory of open quantum systems [1] are found in al- provided by the projection operator techniques [15, 16]
n mostallareasofphysics,rangingfromquantumoptics[2] which are extensively used in nonequilibrium thermody-
a tocondensedmatterphysics[3]andchemicalphysics[4], namics and statistical mechanics [17, 18]. The key con-
u
q from quantum information [5] to spintronics [6]. More- cept of these techniques consists in the introduction of
: over, the theory of open quantum systems provides the a certain projection superoperator P which acts on the
v foundations of quantum measurement theory [7], deco- operatorsof the state space of the total system. The su-
i
X herence [8] and the emergence of thermodynamic behav- peroperator P formalizes the idea of the elimination of
r ior [9]. degrees of freedom from the complete description of the
a states of the total system. Thus, if ρ is the density ma-
In a microscopic approach one regards the total sys-
trix of the composite system, the projection Pρ serves
tem, which is composed of the open system S and its
to represent a simplified effective description through a
environment E, as a closed quantum system following
reduced set of variables. For this reason, the projection
a Hamiltonian time evolution. One of the central goals
Pρ is called the relevant part of the total density ma-
of the theoretical treatment is then the analysis of the
trix,whilethe complementaryprojectionQρ=ρ−Pρis
dynamical behavior of the populations and coherences
referred to as irrelevant part.
which are given by the elements of the reduced density
matrix ρ (t) = tr ρ(t). Here, ρ(t) denotes the density With the help of the projection operator techniques
S E
matrixofthecompositesystemandtr thepartialtrace one derives closed equations of motion for the relevant
E
taken over the environment. part Pρ(t) from which one obtains approximate master
equations by means of a systematic perturbation expan-
IntheMarkovianregimeacompletemathematicalthe-
sion with respect to the system-environment coupling.
ory is available which is based on the concepts of com-
Weshallconcentrateinthispaperonaspecificvariantof
thetechniquewhichisknownasthetime-convolutionless
(TCL) projection operator method [19, 20, 21, 22, 23].
The advantage of this formulation is that it leads to dy-
∗Electronicaddress: [email protected]
†Electronicaddress: [email protected] namicequationsforPρ(t)whicharelocalintimeandin-
‡Electronicaddress: [email protected] volve an explicitly time-dependent generator. A general
2
account of projection operator methods and, in partic- tionsuperoperatorandconnectitwiththemethodofthe
ular, of the TCL approach and its applications may be TCLtechnique. Thisconnectionenablesustodetermine
found in Ref. [1]. higher order corrections in a systematic way and, hence,
Inthestandardapproachtothedynamicsofopensys- to assess the quality of the approximations obtained.
temsonechoosesaprojectionsuperoperatorwhichisde-
fined by the expressionPρ=ρ ⊗ρ , where ρ is some Recently, an entirely different approach has been sug-
S E E
fixed environmental state. This superoperator projects gested,the Hilbertspaceaveragemethod(HAM) [9,10].
the total state ρ onto an un-correlated tensor product This method employs the concept of a best guess for
state. Since ρ is consideredas fixed, it implies that the conditional quantum expectation values. It is based on
E
elements of the reduced density matrix ρ (t) represent the determination of a conditional Hilbert space aver-
S
the relevant variables used for an effective description of age. The methodprovidesuswithasystematicprinciple
the reduced system dynamics. This ansatz for the pro- to estimate quantum expectation values conditioned on
jectionsuperoperatorP iswidely usedinstudies ofopen prescribed values for the expectations of a given set of
quantum systems. It has been used to derive Markovian operators on the total state space. HAM can be used to
andnon-Markovianquantummaster equationsfor many construct effective equations of motion for the given set
applications (see, e. g., [24, 25, 26, 27, 28]). Moreover, of operators and, hence, yields an alternative approach
non-Markovian generalized master equations for the re- to non-Markovian quantum dynamics. It will be shown
duced density matrix have been developed on the basis herethatthemethodoftheHilbertspaceaverageandthe
of phenomenologicalconsiderations [29]. projection operator techniques which are based on the
The paradigm of these and many other approaches is class of correlated projection superoperators are closely
theusageofthereduceddensitymatrixρ (t)asdynami- related. In fact, we are going to demonstrate that HAM
S
calvariableforwhichappropriate(exactorapproximate) represents the lowest order of the TCL expansion corre-
dynamic equations are to be developed. However, it is sponding to this class of superoperators.
important to realize that the projection operator tech-
niquesaremuchmoregeneralandflexible,andthatthey Theapplicationandtheefficiencyofourapproachwill
offer many further possibilities for the construction of be illustrated and discussed here by means of a specific
suitableprojectionsuperoperators. Theonlyformalcon- system-reservoir model. The model consists of a spin
ditionwhichmustbesatisfiedinordertoapplythetech- which interacts with two environmental energy bands
niques is that P is a map which operates on the total through a set of random couplings [40, 41]. This model
state space and has the property of a projection opera- exhibits an unexpected feature. Namely, it turns out
tor,i. e.,P2 =P. This is a verygeneralconditionwhich that the usual Born-Markov approximation fails for this
can be fulfilled in many different ways. model, although the standard Markov condition is satis-
We mention two examples: In the analysis of classi- fied, i. e., although the width of the environmental two-
cal stochastic processes one considers a projectionof the point correlation function is small compared to the re-
form PX(t) = hX(t)i which takes any stochastic pro- laxation time. By contrast, it will be demonstrated by
cess X(t) to its average hX(t)i [19]. With this choice means of a comparison with the numerical solution of
the TCL technique leads to a cumulant expansionof the the full Schr¨odingerequationofthe model thatthe TCL
dynamic equation for the average [30, 31, 32, 33]. In expansion using the correlated projection superoperator
nonequilibrium thermodynamics a further projection of yieldsaccurateresultsalreadyinlowestorderofthe per-
the form Pρ=ρ is often introduced which maps any turbation expansion.
diag
density matrix ρ to its diagonal part ρ in a suitably
diag
chosen basis, yielding the famous Pauli master equation The paper is organized as follows. The class of cor-
in lowest order of the TCL expansion [18]. related projections will be introduced in Sec. II. This
In the present paper we shall construct a class of section also contains a brief general account of the pro-
projection superoperators P which enable the non- jectionoperatortechniques,andintroducesthe basicdy-
perturbative treatment of highly non-Markovian pro- namic equations as well as the perturbation expansion
cesses in open quantum systems. These superoperators of the master equations for the relevant variables. The
project the state of the total system onto a correlated principles and equations of the Hilbert space averaging
system-environment state, i. e., onto a state which con- method are outlined in Sec. III, where we will also de-
tains statistical correlations between certain system and scribetheconnectionbetweentheprinciplesofHAMand
environmentstates. Thus,wegiveuptheparadigmofus- thestructureofthecorrelatedprojectionsuperoperators.
ingthereduceddensitymatrixasthedynamicalvariable, Sec. IV contains the application of the general concepts
and enlarge the set of relevant variables to account for developed here to a specific system-reservoir model. We
statistical correlations which are responsible for strong shalldiscussindetailtheoriginofthefailureoftheBorn-
non-Markovian effects. The idea of introducing addi- Markovapproximation,and show that and why the new
tional variables has been realized in different ways and projectionsuperoperators yield an efficient and accurate
usedinvariouscontexts[34, 35,36,37, 38,39]. Here,we approximation of the dynamics. Finally, we draw our
implementthis ideadirectlyinthe definitionofaprojec- conclusions in Sec. V.
3
II. PROJECTION OPERATOR TECHNIQUES The first equation is just the definition of the reduced
density matrix which is obtained by taking the partial
A. Projection superoperators trace over the environment. The second equation states
that, inorderto determineρ , wedo notreallyneedthe
S
full density matrix of the total system, but only its pro-
We consider an open quantum system S with state
jection Pρ. Hence, the reduced density matrix ρ (t) is
space H which is coupled to an environment E with S
S
foundbytakingtheenvironmentaltraceoftheequations
state space H . The Hilbert space of states of the
E
of motion for Pρ(t).
composite system is given by the tensor product H =
Withinthestandardapproachtothedynamicsofopen
H ⊗H . We assume that the dynamics of the total
S E
systemsusing projectionoperatortechniques onedefines
density matrix ρ(t) of the composite system is governed
a projection superoperator of the form:
bysomeHamiltonianoftheformH =H +V,whereH
0 0
generatesthefreetimeevolutionofthesystemandofthe
Pρ=(tr ρ)⊗ρ , (2.5)
environment, and V describes the system-environment E E
coupling. We work in the interaction representation and
where ρ is some fixed environmental density matrix
writethevonNeumannequationofthecombinedsystem E
referred to as the reference state. This superoperator
as
clearly satisfies our basic conditions (2.3) and (2.4). By
d
ρ(t)=−i[V(t),ρ(t)]≡L(t)ρ(t), (2.1) use of this P the total state of the system is represented
dt
by means of the tensor product state ρ ⊗ρ . Regard-
S E
whereV(t)isthe Hamiltonianinthe interactionpicture, ing the reference state ρ as fixed, one uses the reduced
E
and L(t) denotes the corresponding Liouville superoper- densitymatrixρ (t)asthedynamicalvariable. Applying
S
ator. theprojectionoperatortechniqueonethenfindsamaster
The projection operator techniques are based on the equation for the reduced density matrix ρ whose coef-
S
introduction of a projection superoperator P. This is a ficients are given by certain multitime correlation func-
linear map tions defined by averages with respect to the reference
state ρ . In particular, the master equation obtained
ρ7→Pρ (2.2) E
in second order of the coupling yields the Born approxi-
whichtakes anyoperatorρ onthe totalstate spaceH to mation of the dynamics which involves certain two-time
an operator Pρ on H, and which has the property of a environmental correlation functions.
projection operator: Our new class of correlated projection superoperators
is obtained as follows. We take any orthogonal decom-
P2 =P. (2.3)
positionofthe unitoperatorI onthe statespaceofthe
E
Given a map P with this property one employs the pro- environment,i.e.,acollectionofprojectionoperatorsΠ
a
jection operator techniques to derive from the von Neu- on H which satisfy
E
mann equation (2.1) for the total density matrix ρ(t)
exact and closed equations of motion for its projection Π Π =δ Π , Π =I . (2.6)
a b ab b a E
Pρ(t) (see Sec. IIB). The basic idea underlying this ap- Xa
proach is the following. With an appropriate choice for
the projection superoperatorone intends to obtain a de- Then we can define a linear map by means of
scriptionofthe dynamics ofsystemstateswhichis much
simpler and much more efficient than the description by 1
Pρ= tr {Π ρ}⊗ Π , (2.7)
E a a
means of the full density matrix ρ(t). Thus, the map P N
Xa a
expresses the transition from the full representation in
terms of the total density matrix ρ to a simplified, ef- where N =tr {Π }. It is againeasy to verifythat this
a E a
fective description through a reduced set of dynamical superoperator fulfills the requirements (2.3) and (2.4).
variables defined by the structure of the projection Pρ. By contrast to the standard projection (2.5) which uses
Equation (2.3) is our first basic condition. It is this a representationby a tensor product state, the new pro-
conditionwhichallowsthe formalapplicationofthe pro- jection (2.7) employs a set of un-normalizeddensity ma-
jection operator techniques to open quantum systems. trices
The ultimate goal is to determine form the equations
of motion for Pρ(t) the dynamics of the density matrix ρ(a) =tr {Π ρ} (2.8)
S E a
ρ (t) of the reduced open quantum system. To this end,
S
we need one further condition. Namely, whatever the in order to describe the states of the composite system.
form of the projectionsuperoperatoris, we demand that The set of the density matrices ρ(a)(t) therefore repre-
Pρ contains the complete information requiredto recon- S
sents the dynamical variables defined by the projection
struct ρ . We therefore impose the second basic condi-
S superoperator (2.7). Applying the projection operator
tion:
technique one is then led to a coupled system of equa-
ρ ≡tr ρ=tr {Pρ}. (2.4) tions of motion for the ρ(a)(t), from which one obtains
S E E S
4
thereduceddensitymatrixitselfbymeansoftherelation however, Eq. (2.11) merely implies that the initial state
(2.4): is of the correlated form given by the structure of the
projection superoperator (2.7).
ρS(t)=trE{Pρ(t)}= ρ(Sa)(t). (2.9) In general, the TCL generator K(t) and the inhomo-
Xa geneity I(t) are extremely complicated objects, and the
exact solution of Eq. (2.10) is as difficult as the solution
In the theory of entanglement (see, e. g., Ref. [42])
of the full von Neumann equation for the total system.
a state of the form given by Eq. (2.7) is called separa-
However, Eq. (2.10) can be used as a starting point of
ble or classically correlated [43]. The approachbased on
a systematic perturbation expansion with respect to the
a projection of this form thus tries to approximate the
strengthoftheinteractionHamiltonianV. Withthehelp
total system’s state through a classically correlated but
of the TCL technique one derives a closed expressionfor
non-factorizing state. Of course, one can also construct
the corresponding expansion of the TCL generator:
projection superoperators which lead to non-separable
(entangled) states Pρ. The examples discussed below
∞
belong to the classes of projection operators defined by
K(t)= K (t). (2.12)
n
Eqs. (2.5) and (2.7).
nX=1
The n-th-order contribution is given by
B. Equations of motion
t t1 tn−2
Basically,therearetwovariantsoftheprojectionoper- Kn(t) = dt1 dt2... dtn−1
Z Z Z
atortechnique. ThefirstoneistheprominentNakajima- 0 0 0
×hL(t)L(t )L(t )...L(t )i . (2.13)
Zwanzigmethod[15,16]. Itleadstoafirst-orderintegro- 1 2 n−1 oc
differential equation for Pρ(t) which contains a time
The quantities
integration over the past system history involving a
certain memory kernel. The second variant is known
as time-convolutionless (TCL) projection operator tech- hL(t)L(t1)L(t2)...L(tn−1)ioc ≡
nique [19, 20] which yields a time-local equation of mo- (−1)qPL(t)...L(t )PL(t )...L(t )P...P
i j k
tion for Pρ(t). We shall use this second variant of the X
projection operator technique in the present paper. Its
are known as ordered cumulants [30, 31, 32, 33] and are
advantage is that in any order of the coupling one only
defined by the following rules: (1) Write a string of the
hastosolveafirst-orderdifferentialequationwhichislo-
form PL...LP with n factors of L in between two P’s.
cal in time. It should be emphasized, however, that the
(2) Insert an arbitrary number q of factors P between
general considerations developed here may also be ap-
the L’s suchthat atleast one L stands between two suc-
plied to the Nakajima-Zwanzigprojectionoperatortech-
cessive P’s. The resulting expression is multiplied by a
nique.
factor (−1)q and all L’s are furnished with a time ar-
TheTCLprojectionoperatormethodleadstoanequa-
gument: The first time argument is always t. The re-
tion of motion for the projection Pρ(t) which is of the
mainingL’scarryanypermutationofthetimearguments
general form:
t ,t ,...,t withtheonlyrestrictionthatthe timear-
1 2 n−1
d guments in between two successive P’s must be ordered
Pρ(t)=K(t)Pρ(t)+I(t)Qρ(0). (2.10)
dt chronologically. In the above expression we thus have
t ≥ ... ≥ t , t ≥ ... ≥ t , etc. (3) Finally, the ordered
i j k
This is an exact inhomogeneous linear differential equa-
cumulant is obtained by a summation over all possible
tion of first order. Both the TCL generator K(t) of
insertions of P factors and over all allowed distributions
the linear part and the inhomogeneity I(t) are explic-
of the time arguments.
itly time-dependent superoperators.
In many physicalapplications it may be assumed that
The inhomogeneous part of Eq. (2.10) is determined
the relations
by the projection Qρ(0) of the initial state, where Q =
I−P is the projection superoperator complementary to
PL(t)L(t )...L(t )P =0 (2.14)
P, and I denotes the unit map. We observe that the 1 2n
inhomogeneous term vanishes if the initial state satisfies
hold, which means that any string containing an odd
the relation Qρ(0)=0, i. e., if
number of L’s between successive factors of P vanishes.
Pρ(0)=ρ(0). (2.11) Following the above rules one then finds that all odd-
order contributions K (t) vanish, while the second
2n+1
Obviously,this relationsimplifies theequationofmotion and the fourth-order terms take the form:
(2.10). In the case of the standard projection (2.5) it
impliesthatρ(0)isanun-correlatedtensorproductstate, t
K (t)= dt PL(t)L(t )P, (2.15)
i.e. ρ(0)=ρS(0)⊗ρE. Inthecaseoftheprojection(2.7), 2 Z0 1 1
5
and A and denote it as
t t1 t2 A=Jhψ|Aˆ|ψiK . (3.1)
K4(t) = dt1 dt2 dt3 {hψ|Bˆn|ψi=Bn}
Z Z Z
0 0 0
This expression stands for the average of hψ|Aˆ|ψi over
× PL(t)L(t )L(t )L(t )P
(cid:16) 1 2 3 all |ψi that feature hψ|Bˆn|ψi = Bn but are uniformly
−PL(t)L(t1)PL(t2)L(t3)P distributed otherwise. Uniformly distributed means in-
−PL(t)L(t )PL(t )L(t )P variant with respect to all unitary transformations that
2 1 3
leave hψ|Bˆ |ψi=B unchanged. One may rewrite (3.1)
n n
−PL(t)L(t )PL(t )L(t )P . (2.16)
3 1 2 as
(cid:17)
The performance of the formal expansion outlined A=tr{Aˆαˆ} with αˆ ≡J|ψihψ|K . (3.2)
{hψ|Bˆn|ψi=Bn}
above strongly depends, of course, on the choice of the
projection superoperator P. In other words, the quality How is αˆ to be computed? Any unitary transformation
of the approximation obtained by truncating the expan- that leaves hψ|Bˆ |ψi = B invariant has to leave αˆ in-
n n
sion at a given order n crucially depends on the struc- variant, i. e.:
ture of the chosen projection. It is important to note
that the technique yields an expansion of a certain sys- eiGˆαˆe−iGˆ =αˆ with [Gˆ,Bˆ ]=0. (3.3)
n
tem of equations of motion, and not of the reduced sys-
tem’s density matrix itself. Taking different projection This, however,can only be fulfilled if [Gˆ,αˆ]=0 and this
superoperatorsone uses entirely different sets of dynam- leads to the general form
ical variables which obey completely different equations
ofmotion. Hence,changingtheprojectionsuperoperator αˆ = bnBˆn. (3.4)
amounts to changing the set of dynamical variables and Xn
the whole structure of the equations of motion, and to
(In principal there could be addends of the form, e. g.,
a non-perturbative re-organization of the expansion. It
Bˆ Bˆ , etc., but since all Bˆ we are going to consider
may even happen that the solution of the equations of n m n
below together with zero form a group, those addends
motioninagivenorderforoneparticularprojectionrep-
are already contained in the above sum.)
resents the solution to all orders for another projection.
Furthermoreoneofcoursehasthefollowingconditions:
This point will be illustrated in Sec. IV by means of a
specific system-reservoirmodel.
tr{αˆBˆ }=B . (3.5)
m m
III. HILBERT SPACE AVERAGING
One thus obtains
APPROACH TO THE REDUCED DYNAMICS
B = tr{Bˆ Bˆ }b (3.6)
m m n n
A. The Hilbert space average method Xn
from which the b may be determined. Thus, the con-
TheHilbertspaceaveragemethod(HAM)isinessence n
structionofagivenHilbertspaceaverageis definedwith
a technique to produce guesses for the values of quan-
the helpofEqs.(3.6), (3.4), and(3.2). Accordingto this
tities defined as functions of a wavefunction |ψi if |ψi
scheme“bestguesses”forcertainexpectationvalues will
itself is not known in full detail, only some features of
be produced below. The explanation of HAM in full de-
it. In particular it produces a guess for the expectation
tail is beyond the scope of this text and can be found in
value hψ|Aˆ|ψi if the only information about |ψi is a set
[9, 44].
of different expectation values hψ|Bˆ |ψi = B . Such a
n n
statement naturally has to be a guess since there are in
general many different |ψi that are in accord with the B. HAM, projection operators and dynamics
given set of B but produce possibly different values for
n
hψ|Aˆ|ψi. Thequestionnowiswhetherthedistributionof
In the following we explain how HAM can be used
hψ|Aˆ|ψi’sproducedbytherespectivesetof|ψi’sisbroad
to produce the reduced dynamics of a quantum system
orwhetheralmostallthose|ψi’syieldhψ|Aˆ|ψi’sthatare coupled to some environment, just like the techniques
approximately equal. It turns out that if the spectral described in Sect. II. Consider the full system’s pure
width of Aˆ is not too large and Aˆ is high-dimensional stateatsometimet,|ψ(t)i. LetDˆ(τ)beatimeevolution
almostallindividual |ψiyieldanexpectationvalue close operatordescribingtheevolutionofthesystemforashort
to the mean of the distribution of hψ|Aˆ|ψi’s. In spite of time, i. e., |ψ(t+τ)i = Dˆ(τ)|ψ(t)i. This allows for the
this being crucialfor the following we refer the reader to computation of a set of observables Bˆ at time t+τ:
n
[9]fordetails. Tofindthatmeanonehastoaveragewith
respectto the |ψi’s. We callthis a Hilbert space average B (t+τ)=hψ(t)|Dˆ†(τ)Bˆ Dˆ(τ)|ψ(t)i. (3.7)
n n
6
Now assume that rather than |ψ(t)i itself only the set of Working in the interaction picture the dynamics of the
expectationvaluesB (t)=hψ(t)|Bˆ |ψ(t)iisknown. The full system is controlled by the interaction V(t). The
n n
application of HAM produces a guess for the B (t+τ) time evolution is generated by the corresponding Dyson
n
based on the B (t): series. Thus, assuming weak interactions Eq.(3.12) may
n
be evaluated to second order in the interaction strength
B (t+τ)≈Jhφ|Dˆ†(τ)Bˆ Dˆ(τ)|φiK . (3.8) using an appropriately truncated Dyson series for Dˆ(τ).
n n {hφ|Bˆn|φi=Bn(t)} Thisyieldsafterextensivebutratherstraightforwardcal-
culations for the expectation values corresponding to di-
(Notethatherethe|φiappearratherthanthe|ψ(t)ibe-
agonal elements:
cause those are not actually realized states but denote
the set of states over which the Hilbert space average
B (t+τ)= (3.13)
has to be taken.) Iterating this scheme, i. e., taking the iia
Bn(t+τ) for the Bn(t) of the next step allows for the B (t)+ f(ijab,τ) Bjjb(t) − Biia(t) ,
stepwisecomputationofthe evolutionoftheBn’s. Ifthe iia Xjb (cid:18) Nb Na (cid:19)
set of the B is chosen such that it determines the local
n
state of the considered quantum system completely this
and for the expectation values corresponding to off-
technique produces the local reduced dynamics. The re-
diagonal elements:
sult is of course just like HAM itself only a best guess,
but for appropriate systems this guess can be rather ac-
B (t+τ)= (3.14)
ija
curate.
Here,theBˆn’sarechosenspecificallyasoperatorscor- Bija(t)− 1Bija(t) (f(ikab,τ)+f(kjab,τ)),
responding to elements of the reduced density matrix of 2 Na Xkb
the considered system and the occupation probability of
“energy bands” of the environment: where
Bˆ ≡Bˆ ≡|iihj|⊗Π , (3.9) τ τ′
n ija a f(ijab,τ) ≡ 2 dτ′ dτ′′ (3.15)
Z Z
0 0
where|ii,|jiareenergyeigenstatesoftheconsideredsys-
×tr{Π hi|V(τ′′)|jiΠ hj|V(0)|ii}.
temandΠ isasdescribedinSect.IIAaprojector,pro- a b
a
jectingouttheenergyeigenstatesoftheenvironmentbe-
Those f’s are essentially integrals over the same envi-
longing to an interval ∆E around some mean band en-
a ronmental temporal correlationfunctions that appear in
ergy E labelled by the index a. Let hψ|Bˆ |ψi≡ B ,
a ija ija thememorykernelsofstandardprojectionoperatortech-
then one gets for the elements ρ of the reduced density
ij niques. But here they explicitly correspond to tran-
matrix:
sitions between different energy subspaces of the envi-
ronment. (In Eqs. (3.13) and (3.14) we assumed that
ρ = B . (3.10)
ij jia correlation functions vanish unless they refer to correla-
Xa
tions between parts of the interaction that are adjoints
Thus, the above defined set of expectation values deter- of each other as in Eq. (3.15). We furthermore assumed
mines the reduced state of the system completely. The tr{Πahi|V(τ)|ji} = 0. Both conditions are not neces-
setofstates αˆ whichbelong to the Hilbert spaceaverage sarily fulfilled but apply to the concrete model analyzed
definedbytheB [inthesenseofEq.(3.2)]is,following below.) Those correlation functions typically feature
jia
the scheme described in Sect. IIIA, found to be (short) decay times, i. e, integrating them twice yields
functions which increase linear in time after the corre-
B sponding decay time. Thus, for τ larger than the decay
αˆ = jiaBˆ . (3.11)
N ija time Eq. (3.15) may be written as
Xija a
f(ijab,τ)≈N γ(ijab)τ, (3.16)
(This turns out to be the same state one would have b
gotten from minimizing the purity under the subsidiary
where γ(ijab) has to be computed from Eq. (3.15) but
condition set by the given expectation values B .) A
jia typically corresponds to a transition rate as obtained
comparison with the considerations of Sec. IIA reveals
fromFermi’sGoldenRule. Especiallyitwillonlybenon-
that αˆ has exactly the form produced by the applica-
zero for E −E ≈E −E for otherwise the correlation
i j a b
tion of the projection superoperator P [see Eq. (2.7)].
functionoscillatesrapidlybefore itdecaysandhence the
Exploiting Eqs. (3.1), (3.2), and (3.11) one can write a
corresponding integrals vanish.
specific form of Eq. (3.8) for this case:
InsertingEq.(3.16)intoEqs.(3.13)and(3.14)andas-
suming that the decay times of the correlationfunctions
B (t+τ)≈ (3.12)
ija are small compared to the resulting decay times of the
1
tr |mihl|Π Dˆ†(τ)|iihj|Π Dˆ(τ) B (t). system (which are of the order of 1/γ(ijab)), one can
b a lmb
Xlmb Nb n o transform the iteration scheme into a set of differential
7
equations: model is taken to be H =H +V, where
0
δε
ddtBiia =Xjb γ(ijab)(cid:18)Bjjb− NNabBiia(cid:19), (3.17) H0 = ∆Eσz +Xn1 N1n1|n1ihn1|
δε
d 1 + ∆E+ n |n ihn |, (4.1)
dtBija =−2BijaXkb (γ(kiba)+γ(jkba)). (3.18) Xn2 (cid:18) N2 2(cid:19) 2 2
and
Thissetofdifferentialequationsobviouslydeterminesthe
reduced dynamics of the considered system. It produces V =λ c(n1,n2)σ+|n1ihn2|+h.c. (4.2)
an exponential decay to an equilibrium state. Again nX1,n2
those dynamics are only a guess, but as a guess they
Here and in the following the index n labels the levels
are valid for any initial state regardless of whether it 1
of the lower energy band and n the levels of the up-
is pure, correlated, entangled, etc. In contrast to the 2
per band. σ and σ are standard Pauli matrices. The
standard Nakajima-Zwanzig and TCL methods where z +
overall strength of the interaction is parameterized by
initial states generally produce an inhomogeneity [see
the constant λ. The coupling constants c(n ,n ) are in-
Eq. (2.10)] which may be difficult to handle, HAM al- 1 2
dependent and identically distributed complex Gaussian
lows for a direct guess on the typical behavior of the
random variables satisfying:
system. However, a crucial condition for the application
of the above scheme is that the decay times of the corre-
hc(n ,n )i = 0, (4.3)
1 2
lations are short enough such that even for larger times
hc(n ,n )c(n′,n′)i = 0, (4.4)
the evolution is well described by a Dyson series trun- 1 2 1 2
cated at second order. This means that the scheme will hc(n1,n2)c∗(n′1,n′2)i = δn1,n′1δn2,n′2. (4.5)
break down altogether if the interaction is too strong.
Transforming to the interaction picture we get the von
Neumann equation (2.1) with the interaction picture
Hamiltonian:
IV. APPLICATION
V(t)=σ B(t)+σ B†(t), (4.6)
+ −
A. The model where
B(t)=λ c(n ,n )e−iω(n1,n2)t|n ihn | (4.7)
Toillustratethe generalconsiderationsofthe previous 1 2 1 2
sectionsweinvestigatethemodelofatwo-statesystemS nX1,n2
which is coupled to an environment E [9]. The environ- and
mentconsistsofalargenumberofenergylevelsarranged
in two energy bands of the same width δε. The levels of n2 n1
ω(n ,n )=δε − . (4.8)
1 2
each band are equidistant. The lower energy band con- (cid:18)N2 N1(cid:19)
tainsN levels,theupperbandN levels. Thetransition
1 2
of the two-statesystemis in resonancewith the distance
∆E between the bands (see Fig. 1). B. The standard approach
N 1. Projection superoperator
2
δǫ
In the standard approachone uses a projection super-
operator of the form given by Eq. (2.5). Let us denote
∆E ⊗ the projection onto the lower (upper) band by Π1 (Π2):
N
1
V Π = |n ihn |, (4.9)
1 1 1
δǫ Xn1
Π = |n ihn |. (4.10)
2 2 2
·· system environment Xn2
We consider initial states for which only the lower band
FIG. 1: a two-state system coupled to an environment con- isoccupied: ρ(0)=ρS(0)⊗Π1/N1. Hence,ifwetakethe
sisting of two energy bandswith a finitenumberof levels. reference state
1
ρ = Π , (4.11)
The total Schr¨odinger picture Hamiltonian of the E N 1
1
8
we have This yields the second-order master equation for the re-
duced density matrix:
1
Pρ=(tr ρ)⊗ρ =ρ ⊗ Π , (4.12)
E E S 1
N d 1
1
ρ (t)=γ σ ρ (t)σ − {σ σ ,ρ (t)} . (4.19)
S 2 − S + + − S
dt (cid:20) 2 (cid:21)
and
This is a quantum Markovian master equation in Lind-
Pρ(0)=ρ(0). (4.13)
blad form, where the quantity γ represents the Marko-
2
vian relaxation rate.
In the following we write the elements of the reduced
On the ground of the second order approximation one
density matrix as
couldnaivelyexpectthatthemasterequation(4.19)pro-
ρ (t)=hi|ρ (t)|ji, i,j =0,1, (4.14) videsareasonableapproximationofthereducedsystem’s
ij S
dynamics if the relaxation rate γ is small compared to
2
where |0i and |1i denote the lower and the upper level the band width:
of the two-state system, respectively. If follows from
Eq.(4.13)thattheinhomogeneoustermoftheTCLmas- γ2 ≪δε. (4.20)
ter equation (2.10) vanishes. It can also be verified eas-
However, we are going to demonstrate that this is not
ily with the help of the above forms for the projection
true. Acomparisonwithnumericalsimulationsofthefull
superoperator and the interaction Hamiltonian that the
Schr¨odinger equation and with the prediction of HAM
condition (2.14) holds true. Thus, the TCL generator
shows that the long time dynamics is not correctly re-
K(t) contains only the contributions from even orders of
produced by this master equation.
the coupling strength λ.
1
HAM
2. TCL master equation of second order 0.8 TCL4
Schro¨dinger
From the expression (2.15) for the second-order con- 0.6
tribution of the TCL generator we find: t)
(
1 0.4
t ρ1
K2(t)Pρ(t) = Z dt1f2(t,t1) 0.2 TCL2
0
×[2σ ρ (t)σ −{σ σ ,ρ (t)}]⊗ρ ,
− S + + − S E
0
0 1000 2000 3000
where {·,·} denotes the anticommutator and
t
f (t,t ) = tr {B(t)B†(t )ρ }
2 1 E 1 E
≡ (cid:10)γ h(t−t ) (cid:11) (4.15) FIG. 2: comparison of the numerical solution of the
2 1 Schr¨odingerequationwiththeapproximationsgivenbyHAM
[Eq. (4.21)] and by the second and the fourth order of the
is the two-point environmental correlationfunction with
standard TCL expansion [Eqs. (4.23) and (4.31)]. Parame-
ters: N1 =N2 =500, δε=0.5, and λ=5·10−4.
2πλ2N
2
γ = . (4.16)
2 δε The approximationof HAM leads to the following ex-
pression for the population of the upper level:
The angular brackets in Eq. (4.15) denote the average
overthe randomcouplings c(n1,n2) which is determined ρ (t)=ρ (0) γ1 + γ2 e−(γ1+γ2)t , (4.21)
by use of the relations (4.3)-(4.5). The function h(τ) 11 11 (cid:20)γ +γ γ +γ (cid:21)
1 2 1 2
introduced in Eq. (4.15) is then found to be
whereγ isdefinedbyEq.(4.16)andwehaveintroduced
2
δε sin2(δε·τ/2) a further relaxation rate:
h(τ)= , (4.17)
2π (δε·τ/2)2 2πλ2N
1
γ = . (4.22)
1
wherewehaveassumedaconstantfinitedensityofstates δε
for the environmental energy bands. This function ex-
Ontheotherhand,theTCLmasterequation(4.19)gives:
hibits a sharp peak of width δε−1 at τ = 0 and may
be approximated by a delta function for times t which
are large compared to the inverse band width, i. e., for ρ (t)=ρ (0)e−γ2t. (4.23)
11 11
δε·t≫1 we may approximate:
Thus, the TCL master equation predicts an exponential
f (t,t )≈γ δ(t−t ). (4.18) relaxation of the populations to zero, while the solution
2 1 2 1
9
obtained by means of HAM approaches the stationary be:
population
d 1
ρ (t) = Γ(t) σ ρ σ − {σ σ ,ρ } (4.27)
S − S + + − S
γ dt (cid:20) 2 (cid:21)
ρstat =ρ (0) 1 . (4.24)
11 11 γ +γ 1
1 2 + Γ˜(t) σ σ ρ σ σ − {σ σ ,ρ } ,
+ − S + − + − S
(cid:20) 2 (cid:21)
To judge the quality of the various approximations
where
we have performed numerical solutions of the full
Schr¨odinger equation corresponding to the Hamiltonian Γ(t)=γ (1−γ t), Γ˜(t)=γ γ t. (4.28)
2 1 1 2
definedbyEqs.(4.1)and(4.2). Theinitialstatehasbeen
takentobeoftheform|1i⊗|χi,wheretheenvironmental Equation (4.27) is a master equation with time-
state |χi represents a superposition of the states |n i of dependent relaxation rates Γ(t) and Γ˜(t). The second-
1
the lower band with independent Gaussian distributed order contribution to the rate Γ(t) is given by γ , while
2
random amplitudes of zero mean and equal variances. thefourthorderyieldsthe contributionγ ·(γ t). There-
2 1
For certain parameter ranges we find an excellent agree- fore,thefourth-ordertermoftheTCLgeneratorissmall
mentofthe HAM predictionwiththe simulationresults. compared to the second-order term only if the times t
Anexampleis showninFig.2. Note thatforthe param- considered satisfy the additional condition
etersofthisfigurewehaveγ /δε=3·10−3,suchthatthe
2
γ t≪1. (4.29)
standard Markov condition (4.20) is very well satisfied. 1
We conclude that although the standard Markov con- Thus we see that the occurrence of terms proportional
dition (4.20) is fulfilled the Markovian master equation to t is responsible for strong deviations from the Marko-
(4.19) does not yield a good approximation of the dy- vian behavior. These terms are due to the double-peak
namics for intermediate and long times. In particular, structure of the four-point correlation function f . We
4
its prediction for the stationary state is totally wrong. remark that this structure is markedly different to the
The important point to note is that, in order to judge usual situation of the coupling of an open system to a
the quality of the Markov approximation, an analysis of Bosonic field vacuum, for example. In this case f has
4
the contributions of higher orders is indispensable. We only a single peak and, hence, the above phenomenon of
also note that the same problem occurs if one uses the strongdeviationsfromtheMarkoviandynamicsforweak
Nakajima-Zwanzig master equation. couplings does not occur.
The master equation (4.27) yields the coherences:
ρ (t)=ρ (0)e−γ2t/2, (4.30)
01 01
3. The master equation of fourth order and failure of the
Born-Markov approximation and the populations:
ρ (t)=ρ (0)e−γ2t+γ1γ2t2/2. (4.31)
Tounderstandthe failureofthe Born-Markovapprox- 11 11
imation we investigate the fourth order of the TCL ex-
For γ t≪1 we find the expansion:
1,2
pansion. The contribution of fourth order to the TCL
generator is given by Eq. (2.16). One finds that this 1
ρ (t)=ρ (0) 1−γ t+ γ (γ +γ )t2+... ,
contribution is determined by the two-point correlation 11 11 (cid:20) 2 2 2 1 2 (cid:21)
function (4.15) and by the four-point correlation func-
which is seen to coincide with the corresponding short-
tion:
time expansion of the HAM approximation given by
f (t,t ,t ,t )= tr {B(t)B†(t )B(t )B†(t )ρ } . Eq. (4.21). Thus we conclude that the TCL expansion
4 1 2 3 E 1 2 3 E based on the standard projection reproduces the short-
(cid:10) (4(cid:11).25)
timebehaviorpredictedbyHAMwithinthegivenorders.
Theanalysisshowsthatthisfunctionhastwosharppeaks
Correspondingly,theTCLapproximationoffourthorder
of width δε−1 at t = t , t = t and at t = t , t = t ,
1 2 3 1 2 3
clearly improves the approximation for short times, but
and may be approximated, under the conditions of the
leadsto un-physicalresults forlongertimes anddiverges
previous section, by the expression:
in the limit t → ∞, as is illustrated in Fig. 2. We note
that a similar situation occurs for the spin star model
f (t,t ,t ,t ) ≈ γ2δ(t−t )δ(t −t )
4 1 2 3 2 1 2 3 studied in Ref. [45].
+γ1γ2δ(t−t3)δ(t1−t2), (4.26) Summarizing, the fourth order clearly indicates that
the TCL expansion does not converge uniformly in t. It
where γ1,2 are defined by Eqs. (4.22) and (4.16). only provides a short-time expansion of the dynamics.
Thedouble-peakstructureofthefour-pointcorrelation As a resultofthe emergenceofterms whicharegivenby
expressedby Eq.(4.26) has decisive consequences. With powersoft,itisimpossibletoobtainvalidpredictionson
thehelpofthecorrelationfunctionsgivenabovethemas- the long-time dynamics if one truncates the TCL series
ter equation of fourth order in the coupling is found to at any finite order.
10
C. TCL expansion using the correlated projection Combining Eqs. (4.35) and (4.36) and assuming again
superoperator δε·t≫1 we find the equations of motion:
d γ
In view of the above analysis the following question dtρ(S1)(t) = γ1σ+ρ(S2)σ−− 22{σ+σ−,ρ(S1)},(4.37)
arises: Isitpossibletoconstructanew projectionsuper-
d γ
operator P whose corresponding TCL expansion yields ρ(2)(t) = γ σ ρ(1)σ − 1{σ σ ,ρ(2)}.(4.38)
dt S 2 − S + 2 − + S
the full prediction of HAM already in lowest order, and
leads to a systematic expansion around HAM in higher This is a coupled system of first-order differential equa-
orders? Toanswerthisquestionweconsiderthefollowing tionsforthetwodensitymatricesρ(1)(t)andρ(2)(t). The
S S
projection superoperator: elements of these matrices are written as
1 1 ρ(1)(t)=hi|ρ(1)(t)|ji, ρ(2)(t)=hi|ρ(2)(t)|ji. (4.39)
Pρ = tr {Π ρ}⊗ Π +tr {Π ρ}⊗ Π ij S ij S
E 1 1 E 2 2
N N
1 2
1 1
≡ ρ(1)⊗ Π +ρ(2)⊗ Π . (4.32)
S N 1 S N 2 Schr¨odinger
1 2 1
new TCL2
Thisprojectionbelongstotheclassofsuperoperatorsin- 0.9
troduced in Eq. (2.7). By this ansatz the total system’s HAM
0.8
state is approximatedby a separablebut non-factorizing
state. The dynamical variables are the un-normalized 1 0.7
1
density matrices ρ(1) and ρ(2) which are correlated with ρ 0.6
S S
the projections onto the lower and the upper band, re-
0.5
spectively. The reduced density matrix of the two-state
system is found with the help of Eq. (2.9): 0.4
0 100 200 300 400
ρ (t)=tr {Pρ(t)}=ρ(1)(t)+ρ(2)(t). (4.33) t
S E S S
Assuming that the initial state is of the correlatedform FIG. 3: comparison of the second-order TCL approximation
using the new projection superoperator [Eq. (4.51)], of the
1 1 approximation given by HAM [Eq. (4.21)], and of the nu-
ρ(0)=ρ(S1)(0)⊗ N Π1+ρ(S2)(0)⊗ N Π2, (4.34) merical solution of the Schr¨odinger equation. Parameters:
1 2 N1 =N2 =500, δε=0.5, and λ=0.001.
we have Pρ(0) = ρ(0) and the inhomogeneous term of
Theequations(4.37)and(4.38)cannowbeusedtode-
the TCL equation (2.10) vanishes.
riveanequationofmotionforthereduceddensitymatrix,
making use of Eq. (4.33). First, we get from Eq. (4.38):
1. The second-order master equation dρ(2)(t) = 0, (4.40)
dt 11
d
Using the new projection superoperator (4.32) we get ρ(2)(t) = γ ρ(1)(t)−γ ρ(2)(t). (4.41)
the following second-order TCL equation: dt 00 2 11 1 00
We assume again that initially only the lower band is
d 1 1
Pρ(t) = ρ˙(1)(t)⊗ Π +ρ˙(2)(t)⊗ Π populated:
dt S N 1 S N 2
1 2
= K (t)Pρ(t), (4.35) ρ(2)(0)=0. (4.42)
2 S
It thus follows from Eq. (4.40) that
where the TCL generator takes the form:
ρ(2)(t)≡0. (4.43)
K (t)Pρ(t)= (4.36) 11
2
t From Eq. (4.37) we find
dt h(t−t )
1 1
Z d
0 ρ(1)(t)=γ ρ(2)(t)−γ ρ(1)(t). (4.44)
× 2γ σ ρ(2)σ −γ {σ σ ,ρ(1)} ⊗ 1 Π dt 11 1 00 2 11
1 + S − 2 + − S N 1
h i 1 From Eqs. (4.44) and (4.41) we see that the quantity
t
+ dt1h(t−t1) ρ(111)(t)+ρ(020)(t) is constant. With the help of the initial
Z0 condition (4.42) we thus have
1
×h2γ2σ−ρ(S1)σ+−γ1{σ−σ+,ρ(S2)}i⊗ N2Π2. ρ(020)(t)=ρ(111)(0)−ρ(111)(t). (4.45)