Non-Markovian dynamics without using quantum trajectory Chengjun Wu, Yang Li, Mingyi Zhu, and Hong Guo∗ CREAM Group, State Key Laboratory of Advanced Optical Communication Systems and Networks (Peking University) School of Electronics Engineering and Computer Science, Peking University, Beijing 100871, China (Dated: January 25, 2011) Open quantumsystem interacting with structured environment is important and manifests non- Markovianbehavior,whichwasconventionallystudiedusingquantumtrajectorystochasticmethod. In this paper, by dividing the effects of the environment into two parts, we propose a determin- istic method without using quantum trajectory. This method is more efficient and accurate than stochasticmethodinmostMarkovianandnon-Markoviancases. Wealsoextendthismethodtothe generalized Lindblad master equation. 1 1 PACSnumbers: 03.65.Yz,42.50.Lc 0 2 Whenanopenquantumsysteminteractswithenviron- unitary evolution of the states and the probability flow n ment, it experiences decoherence and dissipation which between these states. Moreover, we also extend this a J leadtolossofinformation. Suchopenquantumsystemis approach to the generalized Lindblad master equation depicted by a reduced density matrix which shows non- whichcandealwithsomestrongcouplingcases[19]. The 3 2 unitaryevolution. Ontheotherhand,theenvironmentis algorithmandnumericalefficiencyaregiven,whichshow classified as Markovianwith no memory effect, and non- thatourmethodismoreefficientandaccuratethanthose ] Markovianwithmemoryeffect. InMarkoviancase,since based on stochastic simulation in most Markovian and h there is no memoryeffect, the quantumtrajectorybased non-Markoviancases. p Monte Carlo wave function (MCWF) method [1–3] and The dynamics of the non-Markovian system is gov- - t quantumstatediffusion(QSD)method[4,5]areapplied. erned by the following master equation [8] n a However, in non-Markovian case, due to memory effect, 1 u theinformationlostbythesystemduringtheinteraction ρ˙(t)= [H ,ρ(t)]+ γ (t)C (t)ρ(t)C†(t) q with the environment will come back to the system in a i~ s Xj j j j 2 [ liaotresrthtiamneMaanrdkosvoiashnocwassem. uch more complicated behav- − 21 γj(t){ρ(t),Cj†(t)Cj(t)}, (1) v Non-Markovian systems are important for their ap- Xj 1 plications to many fields of physics, such as quan- where H is the systemHamiltonian including the Lamb 7 s tum information processing [6, 7], quantum optics [8], shift, C (t)arethe jumpoperatorswhichinduce changes 7 j 0 solid state physics [9], and chemical physics [10]. Re- [e.g., jump from state ψα(t) to ψα′(t) i.e., |ψα′(t)i = . cently, non-Markovian behaviors have also been studied Cj(t)|ψα(t)i/||Cj(t)|ψα(t)i|| ] in the system, and γj(t) 2 in biomolecules where the molecules are embedded in are the decay rates which may take negative values for 1 a solvent and/or in a protein environment [11]. Since some time intervals. The reduced density matrix can be 9 0 there is no true pure state quantum trajectory due to written as [17] : the memory effect [12], the quantum trajectory based v Markovianmethods do not work. Thus, doubled Hilbert Neff i ρ(t)= p (t)|ψ (t)ihψ (t)|, (2) X space (DHS) method [13], triple Hilbert space (THS) α α α r method [14], non-Markovian QSD method [15, 16], and αX=1 a non-Markovianquantum jump (NMQJ) method [17, 18] where p (t) is the probability of the system being in the α areproposedtosolvethenon-Markoviandynamicsofthe state |ψ (t)i at time t. Further, it should be pointed α system where the memory effect is taken into account. out that the effective number of the states N is de- eff However,inordertoobtainhighaccuracy,allthesemeth- termined by C (t)’s [18], Neff p (t) = 1 and that the j α=1 α ods, which are based on stochastic simulations, need to state |ψα(t)i is normalizedP. fulfill a large number of realizations and is very time- To solve the dynamics ofthe system, one shouldknow consuming. So, new methods which are more efficient thetimeevolutionof|ψ (t)ianditsprobabilityp (t). In α α and accurate are highly desired. ourmethod,thetimeevolutionofthestate|ψ (t)iisthe α In this paper, a deterministic method without us- same as that in NMQJ [17]. In NMQJ, the probability ing quantum trajectory is proposed to solve the non- p (t) is calculated in a stochastic way by using quan- α Markovian dynamics. The influence of the environment tum trajectory to N ensemble members. In our method, on the system is divided into two parts, i.e., the non- however, the evolution of probability p (t) is given in a α deterministic way: ′ p˙α(t)=− Γjα(t)pα(t)+ Γjα′(t)pα′(t), (3) ∗Correspondenceauthor:[email protected] Xj (Xα′,j) 2 where Γj(t) =γ (t)kC (t)|ψ (t)ik2 and ′represents where Ni is determined in the same way in Eq. (2) α j j α eff (αP′,j) by taking all the jump operators Rij,s and all the states the summation over all the pairs (α′,j) satisfying ψα(t),s in each ρ (t) into consideraνtion. |ψα(t)i = Cj(t)|ψα′(t)i/kCj(t)|ψα′(t)ik. One finds that jThe evolutionjof state |ψα(t)i is governedby the non- the probabilityofthe state, p (t), changesvia the mech- i α linear differential equation[21] anism of jumps for “out” (α → α′) and “in” (α′ → α), respectively. d i |ψα(t)i=Gˆ(ψα)(t)|ψα(t)i, (8) The numerical simulation corresponding to Eq. (3) is dt i i i straightforward: where Gˆ(ψα)(t) = Hi − i Rji†Rji + i 2 ν ν pα(t+δt)=pα(t)−δt Γjα(t)pα(t)+δt ′Γjα′(t)pα′(t). i Rji|ψα(t)i 2. By combining EqPjsν. (6), (7), Xj (Xα′,j) 2 ν i (4) (8P)jνa(cid:13)(cid:13)nd noting th(cid:13)(cid:13)at ψ1(t) ψ1(t) , ψ2(t) ψ2(t) , ···, Note that there is no stochastic noise and no need to i i i i Ni Ni (cid:12) (cid:11)(cid:10) (cid:12) (cid:12) (cid:11)(cid:10) (cid:12) consider the sign of the decay rate during the simula- ψ eff(t) ψ eff(t(cid:12)) are linear(cid:12)ly(cid:12) independent(cid:12), the (cid:12) i (cid:29)(cid:28) i (cid:12) tion. Additionally, the pα(t)’s in our method do repre- (cid:12)(cid:12)evolution of pα(t) is g(cid:12)(cid:12)iven by sent the probability of the system actually being in the (cid:12) i (cid:12) corresponding pure state ensemble. p˙α(t)=− Γji pα(t)+ ′Γij pα′(t), (9) Consider a particular transition: |ψα′(t)i = i να i να′ j C (t)|ψ (t)i/||C (t)|ψ (t)i||, then the corresponding Xjν (jX,ν,α′) j α j α probability change takes the form: where Γij = Rij ψα(t) 2 and ′ represents the να ν j p (t+δt)=p (t)−δtp (t)Γj(t), (cid:13) (cid:12) (cid:11)(cid:13) (j,Pν,α′) α α α α (5) summationove(cid:13)rallt(cid:12)hepair(cid:13)s(j,ν,α′)satisfying|ψiα(t)i= pα′(t+δt)=pα′(t)+δtpα(t)Γjα(t). Rij ψα′(t) Rij ψα′(t) . It can be easily seen that ν j ν j (cid:12) E.(cid:13) (cid:12) E(cid:13) by s(cid:12)etting n=(cid:13)1an(cid:12)dtaking(cid:13)the decayratesγ(t) into the When the decay rate γj(t) is positive or negative, the (cid:12) (cid:13) (cid:12) (cid:13) equation, Eq. (9) degenerates to Eq. (3). probability flow is from |ψα(t)i to |ψα′(t)i or reversed. Example 1: Detuned Jaynes-Cummings model.– ThishasbeenmentionedinRef.[17]. However,itismore Consider a system with a two-level atom in a detuned explicit in our method. From Eq. (5), it is clear that, in damped cavity, which is governed by the time convolu- thenegativedecayregion,theamountofprobabilityflow tionless master equation [8] only depends on the target state and the probability of the system being in the target state. This is similar to i the situation in NMQJ [17], where the jump probability ρ˙(t)=− S(t){σ+σ−,ρ(t)} 2 in the negative decay region is proportional to the num- 1 1 ber of particles in the target state. These indicate that +γ(t){σ ρ(t)σ − σ σ ρ(t)− ρ(t)σ σ }. − + + − + − 2 2 the trajectory of a particle in NMQJ can not be inter- (10) pretedastrue trajectorysince thejump processdepends on the status of other particles in the system. Because The spectral density of the cavity is supposed to be of true pure state quantum trajectories do not exist in the Lorentzian profile, i.e., J(ω) = γ0λ2 , where non-Markoviandynamics [12], it is not necessary to cal- 2π[(ω0−∆−ω)2+λ2] ∆ = ω −ω is the detuning between the cavity mode culate p (t) in a stochastic way. 0 c α and the atom. To second order approximation, the Next, we extend our method to the recently proposed Lamb shift and the decay rate take the form [8] S(t) = generalizedLindbladmasterequationwhichcansolvethe γ0λ∆ {1−e−λt[cos(∆t)+ λ sin(∆t)]},γ(t)= γ0λ2 {1− dynamics of some highly non-Markoviansystems[19], λ2+∆2 ∆ λ2+∆2 e−λt[cos(∆t)− ∆sin(∆t)]}. In this model, there is only λ d ρ =−i[H ,ρ ]+ Rijρ Rij†− 1{Rji†Rji,ρ } , one jump operator C = σ− = |gihe|, which is a lower- dt i i i (cid:18) λ j λ 2 λ λ i (cid:19) ing operator. We assume thatρ(0)=|ψ1(0)ihψ1(0)| and Xjλ choose|ψ (0)i=(4|ei+3|gi)/5. Actingthejumpopera- 1 (6) tor on the state |ψ (0)i, we get |ψ (0)i=|gi. According 1 2 where i,j =1,2,··· ,n, H are any Hermitian operators, i to Eq. (4), at time t+δt, the probabilities become andRij areanysystemoperators. Itshouldbeindicated λ n p (t+δt)=p (t)−δtp (t)Γ1(t), that ρ(t)= ρi(t). p1(t+δt)=p1(t)+δtp1(t)Γ11(t), (11) iP=1 2 2 1 1 The ith density matrix is decomposed as: where Γ1(t)=γ(t)|he|ψ (t)i|2. 1 1 Neiff Inthis example, ρee(t)is proportionaltothe energyof ρi(t)= pαi(t)|ψiα(t)ihψiα(t)|, (7) the system and p2(t) represents the probability for one αX=1 photon being in the environment. Although p1(t) and 3 p (t) canbe solvedanalytically,inorderto illustrateour 2 method, we use Eq. (11) to do the simulation. The pa- rameters are chosen as ∆ = 12λ,γ λ = 4,λδt = 0.005. 0 Figure 1 (a) shows explicitly the reversal of the proba- FIG.2: (coloronline)Atwo-statesystemcoupledtoanenvi- ronment consisting of two energy bands. Comparison of our method(bluelong-dashedcurve)andMonteCarlosimulation 4 (with N =10 trajectories, green dash-dot curve) to analyt- FIG. 1: (color online) Dynamics of detuned Jaynes- ical result (red solid curve). The parameters are δǫ = 0.31, Cummings model. The initial state is |ψ1(0)i = (4|ei + γ1=γ2 =1 and time step δt=0.01. 3|gi)/5andtheparametersare∆=12λ,γλ=4,λδt=0.005. (a) The probabilities for the system in states |ψ1(t)i and |ψ2(t)i. (b) The population of the excited state ρee (ini- tiallyhigherline)andtheabsolutevalueofthecoherenceρeg than Monte Carlo simulation method. (initially lower line) with three methods: analytic (red solid According to Eqs. (4), we only need to calculate N eff curve),ourmethod(bluelong-dashedcurve)andNMQJ(with statesandchangetheprobabilitiesdeterministically. The N =104 particles in the system, green dash-dot curve). timecostisalmostdeterminedbythecalculationofN eff states. However,the evolutionofN statesis indepen- eff bility flow. We can see from Fig. 1 (a) and (b) that dentwitheachother,sowecancalculatethemparallelly. when the probability flow gets reversed, the energy and In addition, if the jump operatorscan be representedby coherenceoftheatomincrease. Theseshowexplicitlythe sparse matrixes, we only need to calculate the evolution memory effect that the reduced system restores the in- of the states appearing in the decomposition of ρ(0) and formation lost earlier. In Fig. 1 (b), the result of NMQJ usethejumpoperatorstoobtainotherstates. Moreover, (with N = 104 particles in the system) is also given, since the sign of the decay rate makes no difference dur- which shows that our method is more accurate. ing the simulation, in non-Markovian case, our method Example 2: Application to generalized Lindblad mas- is as efficient as it behaves in Markoviancase. ter equation.– To illustrate our method for this kind of Similar to our method, the NMQJ method [17, 18] equation, we consider a two-state system coupled to an needs to calculate N states. However, in addition to eff environment consisting of two energy bands, each with that, NMQJ has to consider the sign of the decay rates a finite number of evenly spaced levels. This may be and generate N random numbers (N ≫N ) to decide eff viewed as a spin coupled to a single molecule or a single the jump process at each time step δt. Apparently, our particlequantumdot[20]. Byusingtime-convolutionless method is more efficient than NMQJ in any case. projection operator technique, to the second order, the InMarkoviancase,theMCWF[1]andQSD[4]method generalizedLindbladmasterequationtakestheform[21] need to realize a large number of trajectories for every state appearing in the decomposition of ρ(0). When d t ρ = dt h(t−t )[2γ σ+ρ σ−−γ {σ+σ−,ρ }], the number of these trajectories is larger than N , 1 1 1 1 2 2 1 eff dt Z 0 which is always the case, our method is more efficient d t ρ = dt h(t−t )[2γ σ−ρ σ+−γ {σ−σ+,ρ }], than them. In non-Markovian case, the DHS method 2 1 1 2 1 1 2 dt Z0 [13], THS method [14] and non-MarkovianQSD method (12) [15, 16] all introduce additional cost for computational where γ h(t − t ),(i = 1,2), is the environment corre- i 1 efficiencycomparedtoMCWForQSD.However,innon- lation function with h(t) = δεsin2(δεt/2) where δε is the Markoviancase,ourmethodisasefficientasitbehavesin 2π(δεt/2)2 widthoftheupperandlowerenergybands. Thereduced Markoviancase. Thus,when the number ofthese trajec- density matrix for the system is given by ρ = ρ1 +ρ2. tories is larger than Neff, our method is obviously more We assume that ρ (0) = |eihe| and ρ (0) = 0. The pa- efficient than them, too. 1 2 rameters are chosen as δǫ = 0.31 and γ = γ = 1. In As for the accuracy, since there is no statistical noise 1 2 Fig. 2 we compare the results of our method, analytical in our method and the error caused by finite time step solution and Monte Carlo simulation which is based on δt is the same, compared with all the methods based theunravelingofthemasterequation(withN =104tra- on stochastic simulation, our method is more accurate. jectories) [21]. Apparently, our method is more accurate Actually, our method is the limit case when the number 4 of realizations in the stochastic based methods tends to tages in efficiency and accuracy. Additionally, we ex- infinite. tendedthis approachto the generalizedLindbladmaster Inconclusion,bydividingthe influenceofthe environ- equation , which is useful to solve the dynamics of some ment on the system into two parts, i.e., the non-unitary highly non-Markoviansystems. evolutionofthesestatesandtheprobabilityflowbetween them, we propose a deterministic method to solve the This work is supported by the Key Project of the Na- non-Makovian dynamics. Compared with the method tional Natural Science Foundation of China (Grant No. based on stochastic simulation, our method has advan- 60837004). [1] J.Dalibard,Y.Castin,andK.Mølmer,Phys.Rev.Lett. and D. J. Tannor, J. Chem. Phys. 123, 204111 (2005) 68, 580 (1992); and references therein. [2] H. Carmichael, An Open System Approach to Quan- [11] P.Rebentrost,R.Chakraborty,andA.Aspuru-Guzik,J. tum Optics, Lecture Notes in Physics (Springer-Verlag, Chem. Phys.131, 184102 (2009). Berlin, 1993), Vol. m18. [12] H. M. Wiseman and J. M. Gambetta, Phys. Rev. Lett. [3] M.B.PlenioandP.L.Knight,Rev.Mod.Phys.70,101 101, 140401 (2008). (1998). [13] H.-P.Breuer,B.Kappler,andF.Petruccione,Phys.Rev. [4] N. Gisin and I. C. Percival, J. Phys. A 25, 5677 (1992); A 59, 1633 (1999). 26, 2233 (1993); 26, 2245 (1993). [14] H.-P. Breuer, Phys. Rev.A 70, 012106 (2004). [5] I.Percival,QuantumStateDiffusion(CambridgeUniver- [15] W. T. Strunz, L. Di`osi, and N. Gisin, Phys. Rev. Lett. sity Press, Cambridge, England, 2002). 82, 1801 (1999). [6] M.A.Nielsen andI.L.Chuang,QuantumComputation [16] J. T. Stockburger and H. Grabert, Phys. Rev. Lett. 88, andQuantumInformation(CambridgeUniversityPress, 170407 (2002). Cambridge, England, 2000) [17] J. Piilo, S. Maniscalco, K. H¨ark¨onen, and K.-A. Suomi- [7] Y. Li, J. Zhou, and H. Guo, Phys. Rev. A 79, 012309 nen, Phys.Rev. Lett.100, 180402 (2008). (2009). [18] J. Piilo, S. Maniscalco, K. H¨ark¨onen, and K.-A. Suomi- [8] H.-P. Breuer and F. Petruccione, The Theory of Open nen, Phys.Rev. A 79, 062112 (2009). Quantum Systems (Oxford University Press, Oxford, [19] H.-P. Breuer, Phys. Rev.A 75, 022103 (2007). 2002). [20] J.GemmerandM.Michel,Europhys.Lett.73,1(2006). [9] See, e.g., C. W. Lai, P. Maletinsky, A. Badolato, and [21] M.MoodleyandF.Petruccione,Phys.Rev.A79,042103 A. Imamoglu, Phys. Rev. Lett. 96, 167403 (2006), and (2009). references therein. [10] J. Shao, J. Chem. Phys. 120, 5053 (2004); A. Pomyalov