Multiple steady-states in nonequilibrium quantum systems with electron-phonon interactions

# Multiple steady-states in nonequilibrium quantum systems with electron-phonon interactions

## Abstract

The existence of more than one steady-state in a many-body quantum system driven out-of-equilibrium has been a matter of debate, both in the context of simple impurity models and in the case of inelastic tunneling channels. In this Letter, we combine a reduced density matrix formalism with the multilayer multiconfiguration time-dependent Hartree method to address this problem. This allows to obtain a converged numerical solution of the nonequilibrium dynamics. Considering a generic model for quantum transport through a quantum dot with electron-phonon interaction, we prove that a unique steady-state exists regardless of the initial electronic preparation of the quantum dot consistent with the converged numerical results. However, a bistability can be observed for different initial phonon preparations. The effects of the phonon frequency and strength of the electron-phonon couplings on the relaxation to steady-state and on the emergence of bistability is discussed.

The existence of a unique steady-state in strongly correlated quantum systems out-of-equilibrium has been a subject of great interest and controversy. For the case of the Anderson impurity model, it has been argued using the Bethe ansatz that a single steady-state solution exists (1). However, recent calculations of the nonequilibrium current based on time-dependent density functional theory seem to indicate that at long times the system reaches a dynamical state characterized by correlation-induced current oscillations (2). Similarly, questions regarding hysteresis, bistability and the dependence of the steady-state current on the initial occupation have been raised in the context of inelastic transport through nanoscale quantum dots (3); (4); (5); (6); (7); (8); (9).

Addressing the issue of a unique steady-state is a challenging task for theory, as systems exhibiting bistability involve strong electron-electron or electron-phonon correlations. Under these conditions, an exact solution is unavailable, and one has to resort to approximate methods or to numerical techniques. The former are based on either a mean-field approximation or on a perturbative scheme, where the inclusion of higher order corrections is not always clear or systematic, and thus may lead to questionable results. Numerically brute-force approaches, such as time-dependent numerical renormalization-group techniques (10); (11); (12), iterative (13); (14); (15) or stochastic (16); (17); (18); (19) diagrammatic methods, and wave function based approaches (20), have been very fruitful, but are limited in the range of parameters and timescales that can be studied.

In this Letter, we address the problem of a unique steady-state in a system with electron-phonon couplings. To study the nonequilibrium dynamics, we develop an approach based on a reduced density matrix (RDM) formalism, which requires as input a short-lived memory kernel (21). We show that if a steady-state exists then it must be unique, regardless of the initial electronic preparation. However, a bistability can develop for different initial phonon states. The relaxation to steady-state and the appearance of the bistability depends on the phonon frequency and the strength of the electron-phonon couplings. To illustrate this, we combine the approach with the multilayer multiconfiguration time-dependent Hartree method (ML-MCTDH) (22) to numerically converge the memory kernel at short times until it decays, and infer from it the dynamics of the system at all times and the approach to steady-state.

We consider a generic model for charge transport through a quantum dot with electron-phonon interaction. The model is described by the Hamiltonian , where is the system (quantum dot) Hamiltonian with creation/annihilation fermionic operators / and energy , where represents the noninteracting leads Hamiltonian with fermionic creation/annihilation operators /, and represents the phonon bath with creation/annihilation bosonic operators / for phonon mode with energy . The coupling between the system and the baths is given by where is the coupling between the system and the leads with couplings strength , and is the couplings between the system and the phonon bath, where is the electron-phonon couplings to mode .

The coupling strengths were determined by various spectral functions. The dot-leads coupling terms were determined from a spectral function of the form , where a tight-binding form was employed: with and . is the chemical potential of the left (L) or right (R) lead, respectively. Similarly, the electron-phonon couplings were determined from a phonon spectral function , where is taken to be of Ohmic form. The dimensionless Kondo parameter, , determines the overall strength of the electron-phonon couplings where is the reorganization energy (or polaron shift) and is the characteristic phonon bath frequency.

Following the derivation outlined in Ref. (21) for the Anderson impurity model, the equation of motion for the RDM, , is given by

 iℏ∂∂tσ(t)=LSσ(t)+ϑ(t)−iℏ∫t0dτκ(τ)σ(t−τ) (1)

where is the system’s Liouvillian, is a trace over the baths degrees of freedom (leads and phonon baths) and is the full density matrix. In the above,

 ϑ(t)=TrB{LVe−iℏQLtQρ(0)} (2)

depends on the choice of initial conditions and vanishes for an uncorrelated initial state (which is the case discussed below), i.e. when , where and are the system and baths initial density matrices, respectively, and . The memory kernel, which describes the non-Markovian dependency of the time propagation of the system, is given by

 κ(t)=TrB{LVe−iℏQLtQLρB} (3)

where , and .

To obtain , one requires as input the super-matrix of the memory kernel, which can be expressed in terms of a Volterra equation of the second type, removing the complexity of the projected dynamics of Eq. (3):

 Extra open brace or missing close brace (4)

with

 Φ(t)=TrB{LVe−iℏLtρB}. (5)

Evaluating the super-matrix requires a tedious calculation similar to that carried out in Ref. (21). The diagonal matrix elements are given by , where

 φmn(t)=TrB{ρB⟨n|∑ktkd(t)a†k(t)|m⟩}. (6)

The indices , , , and can take the values or , corresponding to an occupied or an unoccupied dot, respectively. The off-diagonal elements of are given by , where

 ψmn(t)=TrB{ρB⟨n|∑ktka†k(t)−d†(t)∑αMα(b†α(t)+bα(t))|m⟩}. (7)

Since both the operator appearing in the equation for and the full Hamiltonian conserve the total particle number, is nonzero only for , in which case it has a simple physical interpretation as the time derivative of the dot population. can be expressed in terms of the sum of the left and right currents , where:

 IL,Rm(t)=−2eℏI∑k∈L,Rtk⟨m|d(t)a†k(t)|m⟩, (8)

and is the electron charge. In addition, one can show that is nonzero only for . As a result, the populations and coherences of are decoupled (21), and if one is interested in the behavior of the populations alone, only is required. The resulting equations of motion for the diagonal elements of are given by:

 ∂∂tσii(t)=−1ℏ2∑m∫t0dτκii,mm(τ)σmm(t−τ). (9)

If has a steady-state solution as , then , and can be replaced by a linear set of equations given by:

 ∑mKimσm=0, (10)

where and . The matrix in Eq. (10) has two eigenvalues. The first is corresponding to and . For a physical steady-state solution, both and must share the same sign. Otherwise, the diagonal elements of cannot both be positive. The other eigenvalue is corresponding to , which can never represent a physical solution. Furthermore, since the steady-state depends only on the value of and and since both are independent of the initial dot population, the steady-state is independent of the initial preparation of the dot occupation and therefore, is unique.

The steady-state solution is unique with respect to the electronic initial preparation. This is a strong statement by itself, but it does not rule out bistability for different initial phonon preparations. To address this question, we combined the formalism described above for the RDM with the ML-MCTDH approach in second quantized representation (SQR) (23); (20). The ML-MCTDH-SQR provides a tool to compute the currents in Eq.(8) numerically exactly. The kernel is then obtained by numerically solving Eq. (4). In comparison, for most model parameters studied in this work, it is practically impossible to obtain converged values for the RDM directly from the ML-MCTDH-SQR, since the time to reach a steady-state solution is significantly longer than the maximum simulation time reachable by the ML-MCTDH-SQR. However, since the memory kernel decays on much shorter timescales compared to the RDM itself (21), it is rather straightforward to calculate it using the ML-MCTDH-SQR and then solve Eq. (9) for the RDM.

To characterize the population dynamics, we start with a factorized initial condition of the form , where determines whether the electronic level is initially occupied/unoccupied, represents the initial density matrix of the phonon bath. Hereby two values of the parameters are considered: (corresponding to a phonon initial state equilibrated with an unoccupied dot) and (corresponding to phonons equilibrated to an occupied dot). determines the initial density matrix for the leads. In the above is the inverse temperature. In all results shown below we take and eV.

In Fig. 1 we plot the time evolution of (lower panels) and the corresponding nonzero elements of the memory kernel (upper panels), for two different initial vibrational preparations. We show the time evolution of for different values of the cutoff time at which we assume that the memory kernel has essentially decayed to zero, such that it can be safely truncated. For , it is safe to truncate the memory kernel at while requires a larger cutoff time of . Comparing the time it takes for the memory kernel to decay (upper panels of Fig. 1) with the time taken by the RDM to reach steady-state (corresponding lower panels of Fig. 1), it is clear that the latter is slower by nearly an order of magnitude and in some cases even more. Since the calculation of the memory kernel using the ML-MCTDH-SQR is by far the most time consuming portion of the calculation, the combination with the RDM formalism provides a significant saving, and more importantly extends the ML-MCTDH-SQR approach to regimes inaccessible by direct application.

The inset of Fig. 1 shows the steady-state value of as a function of for . For large values of (short cutoff times) we find that the formalism may lead to unphysical situations in which becomes negative. Of course, this is expected, since only when the memory kernel has decayed to zero does the cutoff approximation provide a physical meaningful solution. As decreases converges and approaches a plateau value. In the present case of parameters, the steady-state value of computed for the two initial vibrational states roughly coincides. However, the dynamics and timescales to relax to the steady-state are clearly sensitive to the initial vibrational preparation.

In Fig. 2 we plot for four different values of and compare the time dependence for four different initial conditions, corresponding to an initial empty () or occupied () dot and to or . For a given , we find that has a unique steady-state solution regardless of the initial dot occupation. This numerically converged result is consistent with the analytical proof given above. In contrast, for two different initial vibrational states, a clear bistability is observed and the RDM decays to a different steady-state solution depending on the value of . We note in passing that a similar phenomenon is known to occur at equilibrium, for example at in the spin-boson model (24); (25).

The appearance of a bistability is consistent with predictions based on a mean field treatment, which is accurate in the adiabatic limit where  (4). For all four frequencies studied, the adiabatic effective potentials (26) shows two distinct minima (upper panel of Fig. 2), corresponding to two possible stable configurations (27). The height of the barrier between the two minima is independent of the phonon frequency, however, as clearly evident in the figure, the width of the barrier increases as decreases. This implies that the tunneling time between the two configurations also increases as decreases, and thus, the bistability (given by the difference between at steady-state for and ) would increase with decreasing , as indeed is the case.

The transient dynamics and the approach to steady state depends sensitively on the preparation, in particular, whether the initial state is close to equilibrium ( and ) or far from equilibrium ( and , ). For small we observe a rapid decay to steady state when the initial preparation is close to equilibrium, while for the nonequilibrium preparation, the population of the dot decays to steady state on a time scale , where denotes the maximum of leads spectral function. As increases, the dynamics becomes more complex. In particular, a longer time scale in the relaxation of developes, consistent with the appearance of a slow tunneling channel between the two configurations. However, the existence of this channel does not lead to a unique steady state solution for the populations at a finite value of , away from the adiabatic limit.

In Fig. 3 we show for different values of the reorganization energy (polaron shift) and the dot energy . As before, we compare the time dependence of for four different initial conditions. In the upper panel we show the corresponding adiabatic effective potentials (26) for the four values of . For small values of , the bistability clearly disappears (left panels of Fig 3). This is consistent with the adiabatic effective potentials showing a single minimum for (in fact, a crude estimate based on a mean-field approach suggests that below the bistability vanishes for the current parameters). Comparing the relaxation time for and , we find that the latter is slower, particularly for the case of .

When is further increased to (corresponding to a nearly symmetric case where the polaron shift equals ) the relaxation time stretches even more and the system decays to a different steady-state depending on the value of , again consistent with the appearance of two stable configurations in the corresponding adiabatic effective potential (upper panel of Fig. 3). While the RDM shows a distinct bistability, it is interesting to note that this is not the case for the current through the quantum dot (not shown), which has a unique steady state for the symmetric case ().

In the upper right panel of Fig. 3 we show results for the case where and eV. The effective adiabatic potential for this case clearly shows two distinct minima, however, the barrier is lower than and eV. Comparing the two right panels of Fig. 3, we find that as and are decreased simultaneously, the bistability decreases and the timescale to relax to steady-state also decreases, consistent with the adiabatic tunneling picture discussed above.

In this work, we have combined a RDM formalism with the ML-MCTDH-SQR approach to study the nonequilibrium dynamics of a many-body quantum system with electron-phonon couplings. For a generic model, which is widely used to describe phonon-coupled electron transport in quantum dots and single-molecule junctions, we showed that the system may exhibit pronounced bistability. The analysis reveals that the bistability increases for decreasing phonon frequency and depends sensitively on the electron-phonon coupling. This was rationalized in terms of the timescales of phonon tunneling between the two adiabatic configurations, the phonon relaxation time, and the electron relaxation. Based on the RDM formalism, we proved that the bistability is associated with different initial phonon preparations and not with a different initial dot occupations. The present results are expected to be of relevance for the interpretation of recent experiments on charge transport in molecular junction, which showed hysteretic behavior (28); (29); (30), and may facilitate further experimental studies in the field of nanoelectronics.

The authors would like to thank Michael Galperin, Abe Nitzan, and Avi Schiller for useful discussion. EYW is grateful to The Center for Nanoscience and Nanotechnology at Tel Aviv University for a doctoral fellowship. HW acknowledges the support from the National Science Foundation CHE-1012479. GC is grateful Yad Hanadiv Rothschild Foundation for the award of a Rothschild Fellowship. MT and ER are grateful to the Institute of Advanced Studies at the Hebrew University for the warm hospitality. This work was supported by the FP7 Marie Curie IOF project HJSC, by the DFG (SFB 953 and cluster of excellence EAM), and used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

### References

1. B. Doyon and N. Andrei, Phys. Rev. B 73, 245326 (2006)
2. S. Kurth, G. Stefanucci, E. Khosravi, C. Verdozzi, and E. K. U. Gross, Phys. Rev. Lett. 104, 236801 (2010)
3. M. Galperin, M. A. Ratner, and A. Nitzan, Nano Lett. 5, 125 (2004)
4. M. Galperin, A. Nitzan, and M. A. Ratner, J. of Phys.: Condens. Matter 20, 374107 (2008)
5. A. S. Alexandrov and A. M. Bratkovsky, J. Phys. Condens. Matter 19, 255203 (2007)
6. A. S. Alexandrov and A. M. Bratkovsky, Phys. Rev. B 80, 115321 (2009)
7. A. A. Dzhioev and D. S. Kosov, J. Chem. Phys. 135, 174111 (2011)
8. K. F. Albrecht, H. Wang, L. Mühlbacher, M. Thoss, and A. Komnik, Phys. Rev. B 86, 081412 (2012)
9. A. O. Gogolin and A. Komnik, “Multistable transport regimes and conformational changes in molecular quantum dots,” arXiv 0207513
10. F. B. Anders and A. Schiller, Phys. Rev. Lett. 95, 196801 (2005)
11. P. Schmitteckert, Phys. Rev. B 70, 121302 (2004)
12. S. R. White, Phys. Rev. Lett. 69, 2863 (1992)
13. S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Phys. Rev. B 77, 195316 (2008)
14. J. Eckel, F. Heidrich-Meisner, S. G. Jakobs, M. Thorwart, M. Pletyukhov, and R. Egger, New J. Phys. 12, 043042 (2010)
15. D. Segal, A. J. Millis, and D. R. Reichman, Phys. Rev. B 82, 205323 (2010)
16. L. Mühlbacher and E. Rabani, Phys. Rev. Lett. 100, 176403 (2008)
17. P. Werner, T. Oka, and A. J. Millis, Phys. Rev. B 79, 035320 (2009)
18. M. Schiró and M. Fabrizio, Phys. Rev. B 79, 153302 (2009)
19. P. Werner, T. Oka, M. Eckstein, and A. J. Millis, Phys. Rev. B 81, 035108 (2010)
20. H. Wang, I. Pshenichnyuk, R. Härtle, and M. Thoss, J. Chem. Phys. 135, 244506 (2011)
21. G. Cohen and E. Rabani, Phys. Rev. B 84, 075150 (2011)
22. H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003)
23. H. Wang and M. Thoss, J. Chem. Phys. 131, 024114 (2009)
24. H. Wang and M. Thoss, New J. Phys. 10, 115005 (2008)
25. U. Weiss, Quantum Dissipative Systems (World Scientific, Singapore, 1999)
26. A. A. Dzhioev and D. S. Kosov, J. Chem. Phys. 135, 074701 (2011)
27. To obtain the adiabatic potential the multiphonon problem was represented by a single mode with potential . The adiabatic effective potential is then given by the integral over the mean force where is the dot energy at position , and , where is the lesser Green function for a fixed phonon coordinate obtained by solving the nonequilibrium Green function on the Keldysh contour (see Ref. [26]).
28. C. Li, D. Zhang, X. Liu, S. Han, T. Tang, C. Zhou, W. Fan, J. Koehne, J. Han, M. Meyyappan, A. M. Rawlett, D. W. Price, and J. M. Tour, Appl. Phys. Lett. 82, 645 (2003)
29. E. Lörtscher, J. Ciszek, J. Tour, and H. Riel, Small 2, 973 (2006)
30. P. Liljeroth, J. Repp, and G. Meyer, Science 317, 1203 (2007)
104303