# Real-time path integral approach to nonequilibrium many-body quantum systems

###### Abstract

A real-time path integral Monte Carlo approach is developed to study the dynamics in a many-body quantum system until reaching a nonequilibrium stationary state. The approach is based on augmenting an exact reduced equation for the evolution of the system in the interaction picture which is amenable to an efficient path integral (worldline) Monte Carlo approach. Results obtained for a model of inelastic tunneling spectroscopy reveal the applicability of the approach to a wide range of physically important regimes, including high (classical) and low (quantum) temperatures, and weak (perturbative) and strong electron-phonon couplings.

In recent years there has been considerable interest in the study of quantum mechanical systems at the nanometer scale that are driven out of equilibrium. Experimental breakthroughs on transport in molecular junctions have uncovered fascinating behavior in molecular systems far from equilibrium Reed1997 (); Joachim2000 (). Much attention has been paid to the study of the transport through strongly correlated systems in which electron correlations are dominant and lead to interesting physics such as the nonequilibrium Kondo effect Gordon1998 (); Cronenwett1998 (); Park2002 (); Liang02 () and Coulomb blockade Elhassid00 () in quantum dots/molecules, tunneling in a Luttinger liquid McEuen99 (); Yacoby02 (), or inelastic effects induced by electron-phonon interactions Nitzan07 () as probed by inelastic electron tunneling spectroscopy Ho98 ().

Exact theoretical treatment of such many-body systems is rather sparse and includes only a small class of simplified model problems Jauho1994 (); Schiller95 (); Voit96 (). For a general solution mean field equations based on the many-body nonequilibrium Green’s function (NEGF) approach can be formulated Jauho_book () and solved numerically. However, the mean-field NEGF approach is quite limited since in many cases it is based on a perturbative scheme and the inclusion of higher order corrections is not always clear or systematic. Traditional determinant Monte Carlo (MC) approaches based on standard discretized PIs Sugar81 (); Hirsch88 () have been used but are numerically rather expensive and seem to be restriced to imaginary-time calculations. Thus, the development of a general approach suitable for the treatment of nonequilibrium many-body quantum systems remains a grand challenge.

In this letter we present a novel approach aimed at obtaining exact numerical results for various dynamical quantities such as the current , conductance , dot population , etc. in a many-body quantum system that is driven out of equilibrium. We focus on a well studied model of inelastic tunneling spectroscopy Flensberg03 (); Galperin2004 (); Mitra2004 (); Hod06 (), where a quantum dot is coupled to two fermionic reservoirs representing the left and right leads at chemical potentials of and , respectively, and to a bosonic bath representing the phonon environment. Motivated by the success of real time PIMC techniques developed for molecular chains Muehlbacher01 (), we adopt a similar procedure and formulate an exact real time path integral (PI) representation for the dynamical quantity of interest. To reduce the computational complexity we integrate out the fermionic leads and the bosonic environment and obtain expressions for their corresponding influence functionals Feynman63 (); Chen87 (). We develop an adequate Monte Carlo procedure where we propagate the density matrix from an initial factorized arbitrary condition towards a steady-state.

A nonequilibrium quantum dot in a phonon environment can be described by the Hamiltonian

(1) | |||||

The time-dependent current from the left lead onto the dot can be written as

(2) | |||||

where is the initial factorizing preparation, is the interaction picture time evolution operator with , and .

We present the approach for a quantum dot which is initially empty; the expressions for the case of an initially occupied dot (or a mixed preparation) can be obtained straightforwardly, as well as those for the right current, the conductivity, the dot’s population etc. After expressing the time evolution operators in Eq. (2) by virtue of Dyson series, can be written as an infinite sum over time-ordered integrals,

(3) |

Here, . The trace over the leads degrees of freedom can now be performed exactly (see, e.g., Ref. Orland_Negele ()), yielding

(4) |

where is a matrix with elements for or for . Here, () denotes the leads’ lesser (greater) self energy in the time domain. in Eq. (3) represents a point correlation function of the dot-phonon subsystem along the Kadanoff-Baym contour . It is most conveniently evaluated in the framework of Feynman PIs, which allows to integrate out the phonon degrees of freedom exactly Feynman63 (), yielding

(5) |

where denotes the Feynman-Vernon influence functional Feynman63 () summarizing the influence of the phonons on the dot, and denote the corresponding forward and backward branches of the dot path, referring to the propagations and , respectively. Since lacks any terms which could flip the state of the dot from empty () to occupied () or vice versa, the contour-ordered sequence of time points uniquely defines the paths and : Every refers to one occurrence of the dot-state altering interaction part in the Dyson series. Therefore, the define the kink times of the dot path Prokofev98 (), between which remains constant. In this spirit, Eq. (3) closely resembles the instanton expansion of the partition function in the spin-boson model Weissbook ().

While the expressions (4) and (5) for and allow for a rather compact expression of , they introduce retardation effects which are arbitrarily long ranged in time, making analytical progress rather cumbersome. To allow for a numerical evaluation of Eq. (3), the integrals in Eq. (3) are approximated dividing the time axis into discrete steps,

(6) | |||||

where . While this introduces a systematic error of the order of and establishes an upper bound to the sum over the kink numbers, the corresponding errors can be made arbitrarily small by adjusting to correspondingly large numbers. In addition, systematic improvement can be made by using a more accurate integration scheme. Performing the sums over all kink numbers and the corresponding (discretized) kink times is now equivalent to summing over all possible (discretized) dot paths , such that Eq. (3) can be written as a discretized PI,

(7) |

where denotes the number of kinks of the path . Eq. (7) can now readily be evaluated by means of PI (or worldline) MC Muehlbacher01 (); Prokofev98 ().

We now turn to discuss the application of the proposed approach to the model system described by the Hamiltonian (1). In the present approach the effect of the leads is fully determined by the self energies (cf. Eq. (4)), which are defined in terms of in Fourier space Jauho_book (). is taken to be energy independent (wide band limit) with a soft cutoff at : . In all results presented below, , , and or to converge the results. Similarly, the phonon influence functional is completely specified by the phonon spectral density Feynman63 (), . For a single phonon coupled to a secondary bath via a coupling constant , becomes a Lorentzian: Leggett87 (). We note in passing that the proposed simulation scheme does neither depend on the particular form of nor of .

In Fig. 1 we plot the time dependent current and the time dependent dot population for different values of the model parameters for (quantum regime) and (strong coupling). Lower panels show the left () and right () currents and upper panels show the average current and the dot population . The left/right currents are characterized by damped coherent oscillation with a long time exponential decay to a steady state value. The present MC scheme provides converged results for the time-dependent current despite the notorious dynamical and fermionic sign problem. This can be contributed to electronic dephasing induced by the leads and the bosonic bath, as well as to the rather small size of the determinants in Eq. (4) (i.e. half the number of kinks on the combined forward-backward dot path) which allows for very fast MC sampling. While the steady state current can be extracted from an exponential fit to , the average current typically decays much faster to steady state, such that in most cases the steady state value could be obtained as the corresponding plateau value.

The fact that approaches faster to steady state is consistent with other flux-based methods, Peskin03 () and can be rationalized by looking at dynamical fluctuations under equilibrium conditions. This is depicted for the case where and the system decays from an initial factorized density to equilibrium. The left/right currents show pronounced coherent oscillations with finite values even at while their average value vanishes for all times as appropriate for equilibrium conditions.

To analyze the limitations of the approach we have conducted a systematic study for the case where (not shown), which numerically is the most difficult case since decoherence effects arising from electron-phonon coupling are neglected. Comparison of the numerical results and the corresponding analytical solution reveals that the present PIMC scheme provides converged results for a wide range of gate and bias voltages, and for a wide range of temperatures spanning the classical limit down to the quantum regime. The method is accurate as long as decays exponentially to steady state within the time window of . As can be seen in Fig. 1 the decay of is slower for decreasing values of and . The approach is still limited for small values of these parameters, as well as other regimes where the steady state timescale and/or the decay time of coherent oscillations are similarily stretched (e.g. very low temperatures); however, improvement of the MC scheme should, in principle, overcome these limitations.

In Fig. 2 we plot the steady state current as a function of the bias for different electron-phonon couplings (lower panel), different phonon frequencies (middle panel), and for different couplings to a secondary bath (upper panel) for (quantum regime). We compare the results to an approximate method based on a generalization of the single particle approximation Wingreen89 () to include the leading order term of the Fermi sea Flensberg03 ().

When the coupling between the primary oscillator and the secondary phonon bath is small we observe steps in the voltage dependent current at integer values of (middle panel) Nitzan07 (). As the coupling to the secondary phonon bath increases or for an Ohmic spectral density the steps diminish and eventually disappear (upper panel), signifying the wide range of phonon frequencies contained in the spectral density. The results for the steady state current shown in the lower panel of Fig. 2 span a wide range of electron-phonon couplings from the Landauer inelastic case, through the perturbative regime to the strong coupling limit. As increases the value of the current decreases from the Landauer inelastic single-channel value to lower values. Our approach clearly captures elastic effects at all values of as depicted by the lower values of the current and by the steps at twice the frequency of the primary phonon mode.

Comparing the exact numerical real-time PIMC results to the approximate expansion method we find that agreement is quantitative for small values of the voltage corresponding to the first step of . For larger values we observe significant deviations between the exact numerical results and the approximate method for . These deviations signify the importance of high order effects of the Fermi sea to the transport process, which are naturally included in the real-time PIMC approach.

In conclusion, we have developed a novel real-time PIMC approach to study the dynamics in open quantum systems that are driven out of equilibrium. The approach is based on expressing the time evolution by virtue of Dyson series, before reducing the dynamics of the entire system by integrating out the fermionic/bosonic baths and introducing exact influence functionals. The remaining infinite sum over contour-ordered time integrals is then evaluated by PI (worldline) MC. We have applied the approach to study the time-dependent current in a well-studied model of inelastic tunneling spectroscopy, where a quantum dot is coupled to two fermionic leads and to a bosonic phonon bath. Numerical results indicate that the approach is robust and can be used for a wide range of model parameters spanning the classical to quantum limits, a range of experimentally accessible chemical bias, different phonon frequencies, weak to strong electron-phonon couplings, and a wide range of couplings between the primary phonon mode and a secondary phonon bath. Furthermore, the approach provides real-time information for various quantities, including the current, conductance, electronic population, spectral function, and more.

We believe that our approach is capable of resolving several shortcomings found in currently used approaches. First, it is not based on any perturbative treatment and in principle can provide exact numerical results. Second, it yields a compact expression for the dynamics which can be evaluated by the proposed PI worldline method in a numerically very efficient way. Furthermore, the real-time propagation scheme allows to study transient phenomena and timescales and also to include time dependent fields. In addition, since it is based on a MC procedure, an enhancement of the accuracy can be obtained by improving the sampling scheme, as well as by including already existing schemes to sooth the dynamical and fermionic sign problem. Finally, it can be applied to a general many body problem, as long as a stable PI worldline approach can be derived.

We would like to thank Andrei Komnik, Guy Cohen, and Abe Nitzan for stimulating discussions. This work was supported by the Israel Ministry of Science (grant to ER). LM would like to thank the Minerva Foundation for financial support.

## References

- (1) M. A. Reed et al., Science 278, 252 (1997).
- (2) C. Joachim, J. k. Gimzewski, and A. Aviram, Nature 408, 541 (2000).
- (3) D. Goldhaber-Gordon et al., Nature 391, 156 (1998).
- (4) S. M. Cronenwett, T. H. Oosterkamp, and L. P. Kouwenhoven, Science 281, 540 (1998).
- (5) J. Park et al., Nature 417, 722 (2002).
- (6) W. J. Liang et al., Nature 417, 725 (2002).
- (7) Y. Elhassid, Rev. Mod. Phys. 72, 895 (2000).
- (8) M. Bockrath et al., Nature 397, 598 (1999).
- (9) O. M. Auslaender et al., Sceince 295, 825 (2002).
- (10) M. Galperin, M. A. Ratner, and A. Nitzan, J. Phys.: Condens. Matter 19, 103201 (2007).
- (11) B. C. Stipe, M. A. Rezaei, and W. Ho, Science 280, 1732 (1998).
- (12) A. P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B 50, 5528 (1994).
- (13) A. Schiller and S. Hershfield, Phys. Rev. B 51, 12896 (1995).
- (14) Y. P. Wang and J. Voit, Phys. Rev. Lett. 77, 4934 (1996).
- (15) H. Haug and A. P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors (Springer, Germany, 1996).
- (16) R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
- (17) R. M. Fye and J. E. Hirsch, Phys. Rev. B 38, 433 (1988).
- (18) K. Flensberg, Phys. Rev. B 68, 205323 (2003).
- (19) M. Galperin, M. A. Ratner, and A. Nitzan, J. Chem. Phys. 121, 11965 (2004).
- (20) A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B 69, 245302 (2004).
- (21) O. Hod, R. Baer, and E. Rabani, Phys. Rev. Lett. 97, 266803 (2006).
- (22) L. Mühlbacher, J. Ankerhold, and C. Escher, J. Chem. Phys. 121, 12696 (2004).
- (23) R. P. Feynman and F. L. Vernon Jr., Ann. Phys. 24, 118 (1963).
- (24) Y. C. Chen, J. Stat. Phys. 47, 17 (1987).
- (25) J. W. Negele and H. Orland, Quantum Many-Particle Systems (Perseus Books, Reading, MA, 1988).
- (26) N. V. Prokofe’ev, B. V. Svistunov, and I. S. Tupitsyn, Phys. Lett. A 238, 253 (1998).
- (27) U. Weiss, Quantum Dissipative Systems (World Scientific, Singapore, 1999).
- (28) A. J. Leggett et al., Rev. Mod. Phys. 59, 1 (1987).
- (29) M. Caspary, L. Berman, and U. Peskin, Chem. Phys. Lett. 369, 232 (2003).
- (30) N. S. Wingreen, K. W. Jacobsen, and J. W. Wilkins, Phys. Rev. B 40, 11834 (1989).