Time-dependent versus static quantum transport simulations beyond linear response
To explore whether the density-functional theory non-equilibrium Green’s function formalism (DFT-NEGF) provides a rigorous framework for quantum transport, we carried out time-dependent density functional theory (TDDFT) calculations of the transient current through two realistic molecular devices, a carbon chain and a benzenediol molecule inbetween two aluminum electrodes. The TDDFT simulations for the steady state current exactly reproduce the results of fully self-consistent DFT-NEGF calculations even beyond linear response. In contrast, sizable differences are found with respect to an equilibrium, non-self-consistent treatment which are related here to differences in the Kohn-Sham and fully interacting susceptibility of the device region. Moreover, earlier analytical conjectures on the equivalence of static and time-dependent approaches in the low bias regime are confirmed with high numerical precision.
A first principles description of quantum transport at the atomic scale is usually based on the non-equilibrium Green’s function (NEGF) technique Kadanoff and Baym (1989), which reduces to the traditional Landauer transmission formalism Landauer (1957); *Buttiker1986 for coherent transport. The device Green’s function is constructed in an ad-hoc manner from the Kohn-Sham (KS) single-particle Hamiltonian of ground state Density Functional Theory (DFT) including the effects of the leads through suitable self-energies Cuevas and Scheer (2010). While such methods provide qualitative results to understand the transport behavior, there are cases where the theoretical steady state current differs from experimental results by orders of magnitude Lindsay and Ratner (2007). One reason for the discrepancy is that the evaluated transmission function exhibits resonances at the KS single particle energies, which are in general not physical.
Recently, several approaches based on time-dependent DFT (TDDFT) have been put forward Tomfohr and Sankey (2001); *Baer2004; *Burke2005; *Kurth2005; *Aspuru2009; Zheng et al. (2007), not only to treat genuine dynamical problems like AC transport or transient effects, but also to remedy the problematic situation in the steady-state case. For optical absorption spectra, TDDFT excited state energies improve significantly on single-particle energy differences that are obtained from the static ground state KS equations. Since the screening properties of the device depend strongly on the optical gap, one would expect also an improved description of the current.
In fact, dynamical corrections to the conventional DFT-NEGF picture have been predicted in the regime of linear response Evers et al. (2004); *Sai2005; *Vignale2009. Interestingly, these corrections were shown to vanish for approximate xc functionals which are local both in time and space. Such adiabatic functionals, like the adiabatic local density approximation (ALDA a.k.a. TDLDA), are quite common in many implementations of TDDFT for optical spectra due to their numerical simplicity. As shown by Stefanucci et. al. Stefanucci et al. (2007) and Zheng et. al. Zheng et al. (2007), TDLDA leads to the well-known Landauer-Büttiker formula for the steady state current if it can finally be reached. It is thus argued that TDDFT and DFT-NEGF should result in the exact same steady state current if the adiabatic approximation is adopted for the exchange-correlation functional. We thus encounter a paradox.
A related open question is whether sizable differences between TDDFT and DFT-NEGF may occur for local functionals beyond linear response. The interest is here in the full I-V characteristics of the device including the bias region around molecular resonances, which is often used as fingerprint of the contacted molecule in experiments. This regime is difficult to access with analytical methods but amenable to numerical simulations. Beyond linear response, more than one solution for steady state current is conceivable, and this was observed for model systems Sánchez et al. (2006). In such a case, the transient current may settle into different steady currents depending on how the bias voltage is turned on. Also, it was argued that persistent oscillating currents can exist and a steady state may never be reached. This is because a device may have bound states with vanishing line widths and corresponding infinite life times. The implications of such bound states for open system were investigated for a model system Khosravi et al. (2008), and the current through the model device was found to oscillate for as long as the simulation time. To address these issues in this letter, we use a recently developed TDDFT method for open systems Zheng et al. (2007) and compare results with those from DFT-NEGF calculations.
Methodology.— For a two–terminal setup, the entire system consists of a device region (), and semi-infinite left () and right () electrodes as shown in Fig 1a. The region , where electron-scattering events take place, is thus an open electronic system.
To obtain the reduced electronic dynamics of the device, we focus on the equation of motion (EOM) of , the reduced single-electron density matrix for D,
where denotes the KS Fock matrix of D, and the square bracket on the right–hand side denotes a commutator. is the dissipation term due to the th electrode and based on the Keldysh NEGF formalism Kadanoff and Baym (1989), it can be expressed as,
where and denote the time domain retarded and lesser Green’s function in the device region, and and stand for the advanced and lesser self-energies of lead .
The transient electric current through the interface (the cross section separating region from the th electrode) can be evaluated via Zheng et al. (2007):
If the bias potentials approach a constant value asymptotically (, a steady state may develop. It has been shown that the steady state current can then be recast in the form Stefanucci et al. (2007); Zheng et al. (2007):
where are Fermi distribution functions and describe the lead-molecule coupling.
For non-local functionals, additional XC contributions appear in the lead Fermi functions and give rise to an effective bias which is smaller than the physical one Evers et al. (2004); *Sai2005; *Vignale2009; Stefanucci et al. (2007). In contrast, the current has Landauer form in the ALDA and seems to be equivalent to the static DFT-NEGF result. However, the device Green’s function
is constructed from a Hamiltonian that depends on the steady-state density. The latter is not necessarily the same as the density obtained from a self-consistent DFT-NEGF treatment or even the equilibrium one. To our knowledge, a general adiabatic theorem for open systems is not available, though adiabaticity has been confirmed for model potentials Scheitler and Kleber (1988). Such a theorem would prove the equivalence of TDLDA and LDA-NEGF approaches from Eq. (5). Hence, we provide a numerical comparison of both approaches for some realistic molecular devices in this letter.
As shown in Ref. Zheng et al. (2007), the expression for in Eq. (2) may be considerably simplified in the adiabatic wide–band limit (AWBL). It depends finally only on the device density matrix , the applied bias V(t) and the WBL parameters which characterize the lead self-energies 111Alternatively (but not employed in this study), a hierarchical equation of motion approach Zheng et al. (2010) was recently developed, for which the hierarchy terminates exactly at the second tier without any approximation.. In this way the EOM is closed. At , the entire system is assumed to be in thermal equilibrium and the initial density matrix is easily constructed using equilibrium Green’s function (EGF) techniques. Subsequently, bias potentials are applied to the contacts and uniformly shift the lead energy levels. The device density matrix is then propagated according to Eq. (1), where in each time step the Poisson equation is solved to update the device Hamiltonian.
Results.— We report the calculations of I-V curves for two molecular systems, a carbon chain and a benzenediol molecule sandwiched between aluminum leads in the (001) direction of bulk Al. The structures are shown in Fig. 1(b) and 1(c) respectively. Simulations for an Al wire were also performed with similar findings as reported below. For each of the molecular systems, we include explicitly 18 Al atoms of each contact in the device region to ensure smooth potentials. The minimal Gaussian basis set STO-3G is adopted in the calculations. For the propagation, the fourth-order Runge-Kutta method is employed with a time step of 0.02 fs to integrate Eq. (1) in the time domain. The ALDA is adopted for exchange and correlation. The bias voltage is turned on exponentially at with a time constant , according to .
Fig. 2 depicts the transient currents for the two molecular systems at different values of with V. Although the initial current traces differ, the final steady-state current is the same in all cases. This was not obvious, as the history of the applied bias voltage is included in the formalism, although the memory effect of the exchange-correlation functional is absent in the ALDA. The establishment of the steady state is ultrafast ( 10 fs) and occurs in every case. We do not observe persistent current oscillations Khosravi et al. (2008), as the term in Eq. (2) ensures proper dissipation in the leads.
TDLDA simulations were then performed for several bias values in the linear response, low voltage, regime. The asymptotic values of the transient currents () were extracted in each case and compared to static DFT calculations in the LDA using the same WBL approximation for the leads in the insets of Fig. 3(a) and Fig. 4(a). Two kinds of static DFT results are reported. In the first case, the device Green’s function is determined fully self-consistently using the non-equilibrium density (). In the second case, the equilibrium density is employed () and the Poisson equation is solved only once, which yields a linear potential drop across the molecule. For both studied systems all approaches yield identical results. This provides numerical evidence for the earlier analytical arguments Evers et al. (2004); *Sai2005; *Vignale2009 on the equivalence of static and dynamic local DFT approaches in linear response.
In another set of simulations, the covered bias range was extended to 5 V. Finite bias simulations allow for the coverage of the resonant tunneling regime in which the bias window contains one or more of the molecular levels. Inspection of the transmission function (not shown here) shows that transport occurs predominantly through the unoccupied levels in both systems. For benzenediol, the conducting molecular orbitals (MO) closest to the Fermi energy () are located 0.92 eV above and 2.58 eV below , respectively. The carbon chain exhibits metallic behavior with a partially filled lowest unoccupied MO derived level Larade et al. (2001).
The full I-V characteristics of both systems are summarized in Fig. 3(a) and Fig. 4(a). The TDLDA () and LDA-NEGF () current traces are virtually identical over the whole bias range including the region of resonant transmission. This level of agreement supports not only the validity of the Landauer picture in the context of TDLDA (Eq. (4)), but also shows that the TDLDA and LDA-NEGF densities are in fact identical for realistic devices even for strong bias voltage. Another reason to study the high bias regime is the possible occurrence of multiple solutions in the self-consistent LDA-NEGF treatment. As discussed by Sánchez et al. Sánchez et al. (2006), a dynamical approach would single out the most stable steady-state solution in such a situation. Attempts to find such multiple solutions were unsuccessful for both devices. We also performed simulations with increased molecule-lead distance in order to cover a larger parameter space of coupling matrix elements. The transient currents show oscillations with slower decay to the steady-state in these cases, but again, the TDLDA results match the corresponding LDA-NEGF simulations exactly.
Comparing now the TDLDA or LDA-NEGF results with LDA-EGF values and focussing first on the bias region below 2 V, we find significant differences which are more pronounced in the carbon chain. The LDA-EGF current is initially higher, which can be traced back to a movement of the LUMO resonance to higher energies in the self-consistent approaches. In fact, TDDFT explicitly includes the potential change due to the induced density, similar to the random-phase-approximation. In the context of harmonic perturbations in the optical range, this leads to the well known shift of TDDFT excited states away from simple KS single-particle energy differences. In order to quantify this effect, we computed the optical spectrum of the isolated benzenediol and carbon chain using a linear response TDDFT method Yam et al. (2003). Fig. 3(b) and 4(b) depict the results using the full response and a non-self-consistent treatment employing the ground state density, i.e. the KS response. Although in both systems the KS and TDDFT spectra differ strongly, the deviations are higher in the carbon chain with a shift of roughly 3 eV compared to a shift of 0.5 eV for benzenediol. This finding is in line with the transport results for V. For larger values of the bias, the observed current traces are more difficult to interpret since the transmission profile is altered in a non-linear fashion with respect to the equilibrium situation.
Summary.—In this letter, we performed TDLDA transport simulations based on a time propagation of the KS density matrix and compared these to conventional static DFT results in the Landauer picture. The studied systems comprise the ubiquitous benzenediol molecule with small equilibrium conductance and an atomic chain with metallic conduction properties. In the regime of linear response, the analytically predicted equivalence of static and time-dependent DFT formalisms for local xc potentials could be confirmed with high numerical precision. Beyond linear response, TDLDA reproduces the result of LDA-NEGF calculations and deviates from static LDA only when the equilibrium density is employed in the latter. This parallels the energy shifts seen in going from the KS response to the full coupled response in DFT absorption spectra. Although the transient current depends on the specifics of applied bias voltage, a unique steady state is reached for the same eventual bias voltage.
More advanced xc functionals beyond the ALDA are expected to lead to an improvement in various ways. For example, non-local functionals provide sizable dynamical corrections to the Landauer formula Sai et al. (2005). Moreover, functionals which exhibit a proper derivative discontinuity refine not only the energetical position of molecular levels, but also the molecule-lead coupling Toher et al. (2005); *Ke2007. Despite of these short-comings, TDLDA is still useful in the qualitative description of transport problems which require a dynamical treatment, like the transient behavior of molecular devices or photo-switches Yam et al. (2008).
We thank the German Science Foundation (DFG, SPP 1243), Hong Kong Research Grant Council (HKU700909P, 700808P, 701307P), Hong Kong University Grant Council (AoE/P-04/08) and National Science Foundation of China (NSFC 20828003) for support.
- Kadanoff and Baym (1989) L. P. Kadanoff and G. Baym, Quantum statistical mechanics: Green’s function methods in equilibrium and nonequilibrium problems (Addison-Wesley, 1989).
- Landauer (1957) R. Landauer, IBM J. Res. Dev. 1, 233 (1957).
- Büttiker (1986) M. Büttiker, Physical Review Letters 57, 1761 (1986).
- Cuevas and Scheer (2010) J. C. Cuevas and E. Scheer, Molecular Electronics: An Introduction to Theory and Experiment (World Scientific, 2010).
- Lindsay and Ratner (2007) S. M. Lindsay and M. A. Ratner, Advanced Materials 19, 23 (2007).
- Tomfohr and Sankey (2001) J. K. Tomfohr and O. F. Sankey, Phys. stat. sol. (b) 226, 115 (2001).
- Baer et al. (2004) R. Baer, T. Seideman, S. Ilani, and D. Neuhauser, J. Chem. Phys. 120, 3387 (2004).
- Burke et al. (2005) K. Burke, R. Car, and R. Gebauer, Phys. Rev. Lett. 94, 146803 (2005).
- Kurth et al. (2005) S. Kurth, G. Stefanucci, C. O. Almbladh, A. Rubio, and E. K. U. Gross, Phys. Rev. B 72, 35308 (2005).
- Yuen-Zhou et al. (2009) J. Yuen-Zhou, C. Rodríguez-Rosario, and A. Aspuru-Guzik, Phys. Chem. Chem. Phys. 11, 4509 (2009).
- Zheng et al. (2007) X. Zheng, F. Wang, C. Y. Yam, Y. Mo, and G. H. Chen, Phys. Rev. B 75, 195127 (2007).
- Evers et al. (2004) F. Evers, F. Weigend, and M. Koentopp, Phys. Rev. B 69, 235411 (2004).
- Sai et al. (2005) N. Sai, M. Zwolak, G. Vignale, and M. Di Ventra, Phys. Rev. Lett. 94, 186810 (2005).
- Vignale and Di Ventra (2009) G. Vignale and M. Di Ventra, Phys. Rev. B 79, 14201 (2009).
- Stefanucci et al. (2007) G. Stefanucci, S. Kurth, E. K. U. Gross, and A. Rubio, Theo. and Comp. Chem. 17, 247 (2007).
- Sánchez et al. (2006) C. G. Sánchez, M. Stamenova, S. Sanvito, D. R. Bowler, A. P. Horsfield, and T. N. Todorov, J. Chem. Phys. 124, 214708 (2006).
- Khosravi et al. (2008) E. Khosravi, S. Kurth, G. Stefanucci, and E. K. U. Gross, Appl. Phys. A 93, 355 (2008).
- Scheitler and Kleber (1988) G. Scheitler and M. Kleber, Zeit. f. Phys. D 9, 267 (1988).
- (19) Alternatively (but not employed in this study), a hierarchical equation of motion approach Zheng et al. (2010) was recently developed, for which the hierarchy terminates exactly at the second tier without any approximation.
- Larade et al. (2001) B. Larade, J. Taylor, H. Mehrez, and H. Guo, Phys. Rev. B 64, 75420 (2001).
- Yam et al. (2003) C. Y. Yam, S. Yokojima, and G. H. Chen, Physical Review B 68, 153105 (2003).
- Toher et al. (2005) C. Toher, A. Filippetti, S. Sanvito, and K. Burke, Phys. Rev. Lett. 95, 146402 (2005).
- Ke et al. (2007) S. H. Ke, H. U. Baranger, and W. Yang, J. Chem. Phys. 126, 201102 (2007).
- Yam et al. (2008) C. Y. Yam, Y. Mo, F. Wang, X. Li, G. H. Chen, X. Zheng, Y. Matsuda, J. Tahir-Kheli, and A. G. William III, Nanotech. 19, 495203 (2008).
- Zheng et al. (2010) X. Zheng, G. H. Chen, Y. Mo, S. K. Koo, H. Tian, C. Y. Yam, and Y. J. Yan, J. Chem. Phys. 133, 114101 (2010).