# Exact stochastic simulation of dissipation and non-Markovian effects in open quantum systems

## Abstract

The exact dynamics of a system coupled to an environment can be described by an integro-differential stochastic equation of its reduced density. The influence of the environment is incorporated through a mean-field which is both stochastic and non-local in time and where the standard two-times correlation functions of the environment appear naturally. Since no approximation is made, the presented theory incorporates exactly dissipative and non-Markovian effects. Applications to the spin-boson model coupled to a heat-bath with various coupling regimes and temperature show that the presented stochastic theory can be a valuable tool to simulate exactly the dynamics of open quantum systems. Links with stochastic Schroedinger equation method and possible extensions to ”imaginary time” propagation are discussed.

###### pacs:

03.65.Yz, 05.10.Gg, 05.70.Ln## I Introduction

Numerous concepts in our understanding of quantum mechanics have emerged from the understanding and description of a system coupled to an environment: measurement, decoherence, appearance of classical world, irreversible process and dissipation… All these phenomena which are often encompassed in the ”theory of open quantum systems”, bridge different fields of physics and chemistry Weiss (1999); Joos et al. (2003); Breuer and Petruccione (2002). During the past decade, important advances have been made in the approximate and exact description of system embedded in an environment using stochastic methods. Recently the description of open quantum systems by Stochastic Schroedinger Equation (SSE) has received much attention Plenio and Knight (1998); Breuer and Petruccione (2002); Stockburger and Grabert (2002). Nowadays, Monte-Carlo wave-function techniques are extensively used to treat master equations in the weak coupling and/or Markovian limitDalibard et al. (1992); Dum et al. (1992); Gisin and Percival (1992); Carmichael (1993); Castin and Mo/lmer (1996); Rigo and Gisin (1996); Plenio and Knight (1998); Gardiner and Zoller (2000); Breuer and Petruccione (2002).

Large theoretical efforts are actually devoted to the introduction of non-Markovian effects. Among the most recent approaches, one could mention either deterministic approaches like Projection Operator techniques Breuer et al. (1999, 2001) or stochastic methods like Quantum State Diffusion (QSD) Diósi and Strunz (1996); Diósi et al. (1998); Strunz et al. (1999); Strunz (2005) where non-Markovian effects are accounted for through non-local memory kernels and state vectors evolve according to integro-differential stochastic equations. In some cases, these methods are shown to be exact Stockburger and Grabert (2002); Shao (2004).

Recently, alternative exact methodsBreuer (2004a); Lacroix (2005) have been developed to treat the system+environment problem that avoid evaluation of non-local memory kernel although non-Markovian effects are accounted for exactly. However, up to now, only few applications of these exact techniques exist Breuer et al. (2004); Breuer (2004b, a); Lacroix (2005); Zhou et al. (2005); Breuer and Petruccione (2007). In all cases, accurate description of the short time dynamics is obtained but long time evolution can hardly be described due to the large increase of statistical errors with time. Although, the possibility to simulate exactly the dissipative dynamics of open quantum systems is already an important step, the challenge to describe long-time evolution is highly desirable to make the techniques more powerful and versatile.

## Ii Exact stochastic equation for the reduced system density

### ii.1 Introduction

In the present work, starting from the exact stochastic formulation of ref. Lacroix (2005) and projecting out the effect of the environment, an equation of motion for the reduced system dynamics is obtained where the environment effect is incorporated through a stochastic mean-field which turns out to be non-local in time. Advantages of the new stochastic theory for the description of long-time evolution of open quantum systems are underlined. We consider here a system (S) + environment (E) described by a Hamiltonian

(1) |

where and denote the system and environment Hamiltonians respectively while is responsible for the coupling. Here we assume that the interaction writes

(2) |

where and correspond to functions of two sets of operators of the system and environment, respectively. In particular, this definition includes non-linear couplings. For the sake of simplicity, we assume an initial separable density . As will be discussed below this assumption could be relaxed. The exact evolution of the system is described by the Liouville von-Neumann equation Due to the coupling, the simple separable structure of the initial condition is not preserved in time. It has however been realized recently in several works using either SSE or path integral techniques that the exact density of the total system could be obtained as an average over simple separable densities, i.e. . In its simplest version, the stochastic process takes the form Lacroix (2005)

(3) |

where the Ito convention for stochastic calculations is usedGardiner (1985). and denote Markovian Gaussian stochastic variables with zero means and variances

(4) | |||||

(5) |

The average over stochastic paths described by Eqs. (3) match the exact evolution. Indeed, assuming that at time the total density writes , the average evolution over a small time step is given by

(6) |

Using statistical properties of stochastic variables (Eqs. (4-5)), we obtain

(7) | |||||

(8) |

Therefore, the last term simulates the interaction Hamiltonian exactly and the average evolution of the total density over a time step reads

(9) |

which is nothing but the exact evolution. Here, the exactness of the method is proved assuming that the density is a single separable density. In practice, the total density at time is already an average over separable densities obtained along each stochastic path, i.e. . Since Eq. (9) is valid for any density written as , by summing individual contributions, we deduce that the evolution of the total density obtained by averaging over different paths is given by which is valid at any time and corresponds to the exact system+environment dynamics.

### ii.2 Stochastic mean-field dynamics

Here, a slightly modified version of the stochastic process is used. It incorporates part of the system/environment coupling using a ”mean-field” approximation in the deterministic evolution. Following ref. Lacroix (2005), we consider the coupled stochastic evolutions, called hereafter ”Stochastic Mean-Field” (SMF) :

(10) |

where

(11) |

The SMF version also provides an exact reformulation of the system+environment problem. Indeed, the two terms in Eqs. (7-8) now read

and properly recombine to recover equation (9).

#### Properties of the SFM theory

Eqs. (10) have several advantages compared to the simple version (Eqs. (3)). First, this stochastic process automatically insures along the stochastic path. In addition, the inclusion of a ”mean-field” term in the deterministic part will always reduce the statistical dispersion compared to the simple stochastic process given by Eqs. (3). This reduction could be significant if quantum fluctuations of the coupling operators and remain small along each path Lacroix (2005, 2007). Indeed, at any time, a measure of the statistical fluctuations is given by

(12) | |||||

Starting from the total density associated to a pure state, the evolution of over a small time step is directly connected to the average quantum fluctuations of and , i.e.

(13) |

where, we have assumed implicitly that all second moments except those given in equations (4-5) cancel out. For comparison, the growth of statistical fluctuations associated to the stochastic process without mean-field (Eqs. (3)) reads

(14) |

Eq. (13) illustrates that the number of trajectories required to simulate the system dynamics will depend on the importance of quantum fluctuations of coupling operators along each path. In addition, a comparison between Eqs. (13) and (14) illustrates that the introduction of mean-field will always improve the numerical accuracy.

### ii.3 Reduced system evolution

In Eq. (10), the influence of the environment on the system only enters through . One expects in general to simplify the treatment by directly considering this observable evolution instead of the complete evolution. To express , we introduce the environment evolution operator

(15) |

Defining the new set of stochastic variables and through and , the evolution of can be integrated as

(16) | |||||

Introducing also the new variables and defined as and , the stochastic equation on the reduced density reads

(17) |

with and where the source term takes the exact form

(18) |

Here, while and are defined by:

(19) | |||||

(20) |

where the environment expectation values are taken at time , i.e. . The two coupled equations, Eqs. (17-18) provide an exact reformulation of the system evolution if and verify

(21) |

In the following, we simply assume that the first term in Eq. (18) cancels out. Substituting Eq. (18) into Eq. (17), we finally obtain an integro-differential stochastic equation for the exact system evolution where the environment effect has been incorporated through the two memory functions (19-20). The interesting aspect of replacing Eq.(10) by Eqs. (17-18) is that, in many physical situations, one can generally take advantage of specific commutation/anti-commutation properties of as well as flexibility in the noise to obtain an explicit form of the memory functions.

## Iii Application to system coupled to a heat-bath

The method is illustrated for systems coupled to an environment of harmonic oscillators initially at thermal equilibrium and it shows that the present stochastic theory can be a valuable tool to simulate exactly open quantum systems. We take

(22) |

and Breuer and Petruccione (2002). The statistical properties of stochastic variables and specified above do not uniquely define the Wiener process. A simple prescription is to further assume

(23) |

There are several advantages to this choice. First, stochastic calculus are greatly simplified. For instance, using standard techniques for system coupled to heat-bath Feynman and Vernon (1963); Caldeira and Leggett (1983) and Ito stochastic rules, shows that and depend only on the time difference and identify with the standard correlation functions Breuer and Petruccione (2002) (see appendix A):

(24) | |||||

(25) |

where

(26) |

denotes the spectral density. No approximation are made to obtain above equations, therefore the average over different stochastic paths match the exact evolution of the system, including all non-Markovian effects.

### iii.1 Equivalent Stochastic Schroedinger Equation formulation

Several works, based on the influence functional method Stockburger and Grabert (2002); Shao (2004); Yan et al. (2004) have led to similar stochastic equations for the reduced density. For instance, authors of ref. Shao (2004); Yan et al. (2004) use an evolution of where the second term in Eq. (18) is absent. As demonstrated below, this term is of crucial importance for applications. In ref. Stockburger and Grabert (2002) and in refs. Breuer et al. (2004); Breuer (2004b); Lacroix (2005) a stochastic formulation of the exact system+environment is given in terms of the Stochastic Schroedinger Equation (SSE) technique. Thanks to the additional stochastic rules (23), Eq. (17) has also its SSE counterpart, where system densities are replaced by and wave functions evolve according to

(27) |

where the bath effect is again incorporated through the mean-field kernel.

### iii.2 Application to the spin-boson model

We illustrate the proposed technique to the spin-boson model. This model can be regarded as one of the simplest quantum open system coupled to a heat bath Leggett et al. (1987) which could not be integrated exactly. In addition, it has been often used as a benchmark for theories of open quantum systems Leggett et al. (1987); Egger and Mak (1994); Diósi et al. (1998); Weiss (1999); Breuer et al. (2001); Stockburger and Grabert (2002); Shao (2004). The system and coupling Hamiltonians are respectively chosen as

where the are the standard Pauli matrices. In the spin-boson model, the numerical solution of equation (17) for the system density, is equivalent to solving three non-linear coupled equations for the .

Fig. 1 shows examples of dynamical evolution of obtained using Eq. (17) and averaging over stochastic trajectories both for weak (filled circle) and strong (open circle) coupling. Results are compared with the Hierarchical approach proposed in ref. Yan et al. (2004). This deterministic approach provides an alternative a priori exact formulation of open quantum systems dynamics and was originally motivated by numerical difficulties encountered in the stochastic theory proposed in ref. Shao (2004). Such difficulties do not occur in the present simulation and much less stochastic trajectories seem to be needed to accurately describe the dynamical evolution. Only trajectories have been used to obtained Fig. 1 leading to statistical errors close to zero (for comparison see discussion in Zhou et al. (2005)). The computer time for the two figures was less than an hour for the weak coupling case up to several hours for the strong coupling case on a standard personal computer. The difference in computing time comes from the fact that smaller numerical time step should be used as the coupling strength increases to achieve good numerical accuracy, the main difficulty being to properly evaluate time integrals in Eq. (18). Denoting the time step by , and have been used for weak and strong coupling respectively.

In the weak coupling case, results of our stochastic scheme displayed in Fig. 1 (filled circles) perfectly match the result of ref. Yan et al. (2004) (solid line). In contrary to ref. Zhou et al. (2005), statistical errors remain small even for long time evolution. The difference in numerical accuracy can be assigned to the second term in Eq. (18) which turns out to be crucial for numerical implementation. Stochastic simulations for strong coupling parameters (open circles) slightly differ from the results obtained with the hierarchical approach in ref. Yan et al. (2004). The numerical convergence of the stochastic simulation presented in Fig. 1 has been checked. Therefore, the difference might be due to the fact that the numerical accuracy depends on the truncation scheme used in the hierarchy, even though the method of ref. Yan et al. (2004) is exact.

### iii.3 Comparison with the Time-convolutionless method (TCL)

The possibility to simulate exactly the system dynamics can also serve as a benchmark for other methods. For instance, we compared the exact stochastic scheme with the approximate Time-Convolutionless (TCL) projection operator method of ref. Breuer et al. (1999); Breuer and Petruccione (2002). Figure 2 presents the results of the exact stochastic simulation compared with the TCL2 method applied to the spin-boson model in ref. Breuer et al. (1999). In this figure different cases corresponding to either low or high temperature regime and weak or strong coupling are presented. We see that the best agreement is obtained in the weak coupling and high temperature case. In general, the TCL2 method compares well with the exact simulation if the coupling is rather small. As the coupling increases (lower panels of figure 2), the difference between the TCL technique and the exact method increases. Note that the TCL method seems to systematically underestimate the damping. Note finally that the use of TCL4 equations instead of TCL2 does not improve the comparison.

### iii.4 Discussion of approximate system evolution obtained with real noise and Hermitian system densities in stochastic evolution

Applications presented above use specific constraints on the Markovian process given by Eqs. (23). This prescription greatly simplifies stochastic calculus. For instance, simple exact expressions have been obtained for and when the system is coupled to a heat-bath of harmonic oscillators (Eqs. (24-25)). The main consequence of Eqs. (23) is that and should be complex stochastic variables leading to non Hermitian densities along paths. As illustrated above, such a stochastic process could be used to simulate exactly the reduced density evolution. The main disadvantage of non-Hermitian densities is however that system observables could hardly be interpreted. We discuss here the possibility to perform stochastic evolution of reduced Hermitian densities.

Relaxing the constraints given by Eqs. (23), authorizes to choose and as real stochastic variables, which automatically insures that and remain hermitian. This alternative has however two major drawbacks. First, one cannot anymore have an equivalent SSE picture. Second, while still identifies with Eq. (24), no simple exact expression can be worked out for . However, since this kernel is a functional of , a hierarchy of more and more accurate approximations could be obtained by successive replacements of into Eq. (19) by its integral expression, Eq. (16). In the present work, we concentrate on the simplest case where is replaced by in the time integral of the memory kernel. In this limit, also reduces to Eq. (25). By doing this approximation, the stochastic process is not exact anymore. Fig. 3 presents a comparison of the exact stochastic simulation obtain with complex noise (filled circles) and the approximate case with real noise (open circles). The parameters of the spin-boson model are the same as in Fig. 1.

This figure shows that the approximate scheme with real noise is very close to the exact simulation even in the strong coupling limit. In the latter case, only at very large time, the two simulation start to deviate slightly from one another. For completeness, approximate stochastic simulations obtained for cases presented in figure 2 are compared to the exact scheme in figure 4. We see that except for the weak coupling and low temperature case, the approximate simulation is very close to the exact case. It is worth mentioning that approximate simulation presented here uses the simplest prescription for . Therefore, improved description could a priori be obtained using better approximation of obtained with the method described above. This example is very encouraging and provides a new method to simulate open systems with a stochastic process preserving the hermitian properties of the system density.

## Iv Conclusion

The results obtained with the stochastic theory for the spin-boson model are very encouraging. The theory turns out to be accurate not only to simulate short but also long-time dynamics and does not seem to suffer from numerical instability quoted in ref Zhou et al. (2005). It is worth mentioning that optimization techniques proposed in ref. Lacroix (2005) which are not used here, can be implemented to further reduce the number of stochastic paths. Besides the numerical aspects, the link with classical dissipations dynamics could be easily made, similarly to ref. Stockburger and Grabert (2002).

Our stochastic approach could be generalized to any initial conditions that could be written as a mixing of separable densities, i.e. with where the are complex coefficients. Then, the complete exact dynamics is recovered by both averaging over trajectories starting from each individually and averaging over the initial conditions. The theory is also not restricted to real time evolution. Statistical properties of the system+environment could be studied by considering imaginary time propagation, i.e. . Imaginary time propagations lead naturally to densities written as a mixing of separable densities and could then serve as initial conditions for real time evolution. By combining both imaginary time and real time propagation, general physical problems similar to those depicted in ref. Stockburger and Grabert (2002) can be addressed. The main limitation of the technique is clearly the choice of coupling operator which should give simple memory functions (Eq. (19 - 20)). It is however worth mentioning that most of the coupling operators used in the context of open quantum systems already enter into this category Carmichael (1993); Weiss (1999); Gardiner and Zoller (2000); Breuer and Petruccione (2002). We believe that the stochastic theory presented in this paper can be a valuable tool to access exactly the dynamics of more complex open quantum systems. We presented here specific applications on systems coupled to a heat-bath. The approach can however be applied to various types of environments and couplings which might be of great interest to address dissipation, measurement and/or decoherence problems in quantum systems.

###### Acknowledgements.

The author thank G. Adamian, N. Antonenko, S. Ayik, and B. Yilmaz for valuable discussions and C. Simenel, G. Hupin for his careful reading of the manuscript.## Appendix A Proof of Eq. (18) for a heat-bath of Harmonic oscillators

We give here a proof of Eq. (18) where and identify with Eqs. (24-25). The environment is assumed to be a set of harmonic oscillators, labeled by ”” associated to creation/annihilation (,), i.e.

(28) |

If thermal equilibrium is initially assumed, the environment density writes as a product of densities of each oscillator, where each density can be written as Gaussian operators (see for instance Gardiner and Zoller (2000)) determined by first and second moments of the . The time evolution of the environment is given by Eq. (10) where and where the fluctuating variables verify (according to Eq. (23)):

(29) |

Above prescription and the specific form of induces important simplifications listed below. First, the initial product form of the environment density is preserved along the stochastic paths where each oscillator density verifies and where for each pairs of densities . Second, due to the linear coupling operator, the Gaussian nature of initial densities is also preserved along paths. Therefore, we can equivalently solve the density equation of motion or follow first and second moments of each in time. Here we consider the second strategy.

From the evolution, the equation of motion of each pairs and reads:

(30) | |||||

where we have introduced the notation and . Here, , and denote the second moments of the , operators:

According to the stochastic environment dynamics, we can show that these moments simply evolve as

(31) |

Since we assume that each oscillator is initially at thermal equilibrium, we deduce that second moments are constant in time with , while . Here, we have introduced the standard function Gardiner and Zoller (2000): . Substituting in Eqs. (30-A) and using standard integration techniques Breuer and Petruccione (2002) leads finally to

(32) |

Accordingly, each position operator entering in reads:

Assuming the initial conditions , substituting in the expression of and introducing the spectral density (Eq. (26)), we finally recover equation (18) where and are respectively given by Eqs. (24) and (25).

### References

- U. Weiss, Quantum Dissipative Systems (World Scientific, Singapore, 2nd ed., 1999).
- H. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford, 2002).
- E. Joos, H. Zeh, C. Kiefer, D. Giulini, J. Kupsch, and I.-O. Stamatescu, Decoherence and the Appearance of a Classical World in Quantum Theory (Springer, New York, 2003).
- M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
- J. T. Stockburger and H. Grabert, Phys. Rev. Lett. 88, 170407 (2002).
- J. Dalibard, Y. Castin, and K. Mølmer, Phys. Rev. Lett. 68, 580 (1992).
- R. Dum, P. Zoller, and H. Ritsch, Phys. Rev. A 45, 4879 (1992).
- N. Gisin and I. Percival, J. Math. Phys. A25, 5677 (1992).
- H. Carmichael, An Open Systems Approach to Quantum Optics (Lecture Notes in Physics, Springer-Verlag, Berlin, 1993).
- Y. Castin and K. Mo/lmer, Phys. Rev. A 54, 5275 (1996).
- M. Rigo and N. Gisin, Quantum Semiclass. Opt. 8, 255 (1996).
- W. Gardiner and P. Zoller, Quantum Noise (Springer-Verlag, Berlin-Heidelberg, 2000).
- H.-P. Breuer, B. Kappler, and F. Petruccione, Phys. Rev. A 59, 1633 (1999).
- H.-P. Breuer, B. Kappler, and F. Petruccione, Ann. Phys. 291, 36 (2001).
- L. Diósi and W. T. Strunz, Phys. Lett. A224, 25 (1996).
- L. Diósi, N. Gisin, and W. T. Strunz, Phys. Rev. A 58, 1699 (1998).
- W. T. Strunz, L. Diósi, and N. Gisin, Phys. Rev. Lett. 82, 1801 (1999).
- W. T. Strunz, New Journ. of Phys. 7, 91 (2005).
- J. Shao, J. Chem. Phys. 120, 5053 (2004).
- H.-P. Breuer, Eur.Phys.J. D 29, 105 (2004a).
- D. Lacroix, Phys. Rev. A 72, 013805 (2005).
- H.-P. Breuer, D. Burgarth, and F. Petruccione, Phys. Rev. B 70, 045323 (2004).
- H.-P. Breuer, Physical Review A 69, 022115 (2004b).
- Y. Zhou, Y. Yan, and J. Shao, Europhys. Lett. 3, 334 (2005).
- H.-P. Breuer and F. Petruccione, Phys. Rev. E 76, 016701 (2007).
- W. Gardiner, Handbook of Stochastic Methods (Springer-Verlag, Berlin-Heidelberg, 1985).
- D. Lacroix, Ann. of Phys. 322, 2055 (2007), eprint quant-ph/0606246.
- R. P. Feynman and J. Vernon, F. L., Ann. Phys. 24, 118 (1963).
- A. O. Caldeira and A. J. Leggett, Ann. Phys. 149, 374 (1983).
- Y. Yan, F. Yang, Y. Liu, and J. Shao, Chem. Phys. Lett. 395, 216 (2004).
- A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
- R. Egger and C. H. Mak, Phys. Rev. B 50, 15210 (1994).