Higher harmonics and ac transport from time dependent density functional theory
We report on dynamical quantum transport simulations for realistic molecular devices based on an approximate formulation of time-dependent Density Functional Theory with open boundary conditions. The method allows for the computation of various properties of junctions that are driven by alternating bias voltages. Besides the ac conductance for hexene connected to gold leads via thiol anchoring groups, we also investigate higher harmonics in the current for a benzenedithiol device. Comparison to a classical quasi-static model reveals that quantum effects may become important already for small ac bias and that the full dynamical simulations exhibit a much lower number of higher harmonics. Current rectification is also briefly discussed.
Keywords:Time-dependent Density Functional Theory molecular electronics ac transport
Molecular electronics, the use of single molecules as functional entities in nanoscale electronic circuits, received significant interest over the past years (1). The synergy of experimental achievements and first principles atomistic modeling was especially fruitful in this field. Key to a correct interpretation of conduction measurements in break junctions or scanning tunneling microscopy (STM) is for example the precise knowledge of the molecular geometry within the junction. This information can nowadays be obtained from Density Functional Theory (DFT) in a routine fashion. Besides structural data, DFT may also be used to estimate I-V characteristics in conjunction with the Landauer-Büttiker formalism (2); (3); (4).
With the enormous chemical diversity of even small to medium sized molecules, the initial hope was to choose appropriate compounds for a given desired circuit element and fine-tune their electronic structure by suitable functionalization. While this approach met practical problems (5), several groups are now investigating the possibility of using time-dependent fields to gain additional control over the conduction process (6); (7); (8). As an example, application of visible and UV light induces a reversible ring-opening/ring-closure reaction in some organic molecules that leads to sizable on/off ratios of the current (7); (9); (10). Light emission from molecular junctions may also provide further means to characterize the device (11); (12); (13). In many of these experiments, the frequency of the the applied radiation is small compared to the plasmon frequency of the metallic leads. One can then assume that the effect of the time-dependent field is just a rigid shift of the contact energy levels, leading to a sinusoidal bias potential across the junction. The resulting picture of ac conduction raises a couple of interesting questions in the context of molecular electronics. Do conventional organic molecules behave as classical capacitors or inductors, or does their ac response depend on frequency and chemistry in a non-trivial way? This question was already addressed by Fu and Dudley in the 1990s for the model system of a resonant tunneling diode (14), and recently, several groups presented DFT based calculations of ac conductances for realistic molecular junctions (15); (16); (17); (18); (19).
The time-dependence of the voltage has also other measurable consequences: If the I-V characteristic of a device is non-linear, the Fourier transform of the time-dependent current will feature peaks at multiples of the driving frequency . These higher harmonics were utilized in STM to characterize also non-conducting samples (20); (21). With the currently available instrumentation, harmonics at optical frequencies are clearly difficult to detect. However, the generation of higher harmonics implies current rectification (the change of dc current due to an additional ac bias) and is therefore indirectly accessible. In recent experiments the measured non-linearities and the rectified current were used to estimate the field enhancement in metallic nanoconstrictions subjected to visible light (22); (23). On the theoretical side, the modeling of photocurrents was pioneered by Tien and Gordon (24), while Cuevas and coworkers presented a number of detailed DFT based calculations on this topic in recent years (25); (26); (27). A direct estimation of the amplitude and number of higher harmonics in realistic molecular junctions seems to be missing. In principle, Green’s function approaches in the energy domain can be utilized for this purpose, but require truncations in the number of photons absorbed and emitted in the tunneling process (28); (29); (1).
In this article we apply an approximate method based on time-dependent density functional theory to investigate the ac response of several realistic molecular devices. The density matrix of the device is propagated according to a Liouville-van Neumann equation including the effects of dissipation and decoherence due to the macroscopic leads. The latter are characterized by an embedding self-energy in the wide-band approximation. At each time step the potential in the device region is determined by solving the Poisson equation and used to update the time-dependent DFT Hamiltonian in a selfconsistent fashion. The method allows one to compute the time-dependent current for arbitrary time-dependent bias potentials. After a short description of the numerical scheme in Sec. 2, we investigate the frequency dependent admittance for a typical molecular junction in Sec. 3. The remainder of the paper (Sec. 4) is devoted to an analysis of the higher harmonics in a molecular junction based on benzenediol. The goal is here to contrast the results from a fully quantum-mechanical calculation with classical estimates.
The following will outline the approach leading to the presented results. Apart from the Hamiltonian entering the scheme, we follow the method introduced by Zheng and co-workers (30). Our systems are partitioned into a central device region (D) - which also includes two layers of contact material - and semi-infinite periodic leads on either side of the device (L,R). The leads are identical and the whole device is assumed to be in thermal equilibrium ( K) at a common chemical potential for . In consequence, current only flows due to an applied time-dependent bias of arbitrary shape and magnitude. As already mentioned in the introduction, we employ a propagation scheme based on the Liouville-van Neumann equation for the Kohn-Sham density matrix. According to Ref. (30), a closed equation for the reduced density matrix of the device region may be derived, featuring dissipation terms arising from the interaction with the semi-infinite leads (in atomic units):
Here, is the reduced Kohn-Sham density matrix written in a basis of localized atom-centered basis functions . The have been derived in the framework of quantum dissipation theory using the Keldysh non-equilibrium Green’s function formalism. It is possible to write down a formally exact and explicit expression for the dissipation terms using hierarchical equations of motions (31) applicable to arbitrary non-Markovian dissipative systems. However, we aim for numerical efficiency and use the Markovian approximation to proposed in the original article (30) within the wide band limit (WBL):
where denotes the real part of the self-energy and its imaginary part. The self-energy itself is calculated from the device-lead couplings and the surface
Green’s function of the leads. The terms can also be readily
calculated from the Hamiltonian and the applied bias potential. The
initial density matrix is obtained from equilibrium Green’s function theory
using the same input Hamiltonian and the wide band limit for consistency. The equation of motion is then numerically propagated by a Runge-Kutta method.
In each time step, the time-dependent current can then be calculated through
In the original scheme, the Hamiltonian governing the propagation is set to be from time-dependent density functional theory (TD-DFT (32)) in the adiabatic approximation. However,
the systems we investigate involve more than a hundred heavy atoms in the device region, and time step converged trajectories require a resolution of ten attoseconds or less. At present state, it is
numerically extremely expensive to extract information about the
admittance behaviour for such demanding systems by use of
TD-DFT. Consequently, we employ the above propagation scheme in
conjunction with an approximate form of
TD-DFT, namely time dependent density functional based tight-binding (TD-DFTB) (33); (34); (35) in addition to the WBL already mentioned.
The Hamiltonian matrix then takes the following form:
Eq. 4 results from an expansion of the electron density in the device region around a constant reference density , assumed to be a superposition of the atomic densities . The zeroth order represents a DFT Hamiltonian (with the PBE XC-functional (36) entering in our case) evaluated at the reference density in a two-center approximation. Since the evaluation of this term only requires information about the distance between atoms, it can be built from pre-calculated lookup tables almost instantaneously. The second term accounts for deviations from the reference density. At each time step, the charge density in the device region is computed from the density matrix and the basis function overlap . The solution of the Poisson equation with boundary conditions matching the applied bias then determines the potential shifts . A more detailed description of the method may be found in Ref. (35).
3 Admittance of a Au-Hexen-Au nanojunction
Using the above approximations, we are able to treat junctions of considerable size and can move away from mere test systems. We investigated the admittance behaviour of a Hexene molecule coupled to Gold leads using thiol bridges. The leads consist each of 172 Au atoms which set up a bulk by periodical replication and are used to calculate the self-energy. The experimental lattice constant of approximately 4.08 Å (37) has been used. The device region includes two layers of the described contacts which are aligned along the (111)-surface in transport direction with 86 Gold atoms in total. Pyramids consisting of four atoms have been placed on top of the surfaces with the same interatomic distances as in the bulk. The molecule has been optimized in vacuum using DFT with the PBE functional and a 6-31G* basis set. In order to get a good guess for the binding angle, an S-Au group was previously attached to either side. This approximate geometry was then placed symmetrically between the pyramids without further optimizations (see Fig. 1). In our transport calculations, all atoms are described within a minimal basis set: l=0 for H, up to l=1 for C and S, and up to l=2 for Au.
In the time-dependent simulation, we applied a Lorentzian input signal
with fs, fs, V and used its Fourier transform as well as the Fourier transform of the current response to obtain a complex admittance curve according to . The result for the Au-Hexen-Au junction is depicted in Fig. 2a). Qualitatively, one can identify capacitive behaviour for small frequencies in analogy to a classical RC-circuit (up to second order) (15):
In an earlier study, we observed the quadratic increase of the conductance (Re[Y]) and the linear decrease of the susceptance (Im[Y]) also in a simpler system with little transmission in the vicinity of the Fermi energy (38).
For a single level model, the admittance may be determined analytically as shown by Fu and Dudley (14). Given a dc Breit-Wigner transmission of the form
where is the level energy, the ac admittance reads
where and denotes the quantum of conductance. In Fig. 2b), we plot a visualization of these formulas for different model parameters. If the conduction is off-resonant in the dc limit, a capacitive regime with negative values of Im[Y] occurs at finite , while resonant transmission leads to inductive behaviour. Generally, a broader peak in the transmission also leads to a significantly greater effect on the shape of the admittance, especially its real part. Resonances at energies far from contribute less to the admittance. Consequently, one expects to be able to relate the profile of the admittance to broad features in the dc transmission, at least close to the Fermi energy.
We therefore computed also the transmission for the Au-Hexen-Au junction in the Landauer-Büttiker framework at the DFTB level (4). In Fig. 2c), one observes a large number of peaks with a width close to zero in the energy range [-4,-2] eV. These likely correspond to localized metal states that do not couple to the left and right lead and therefore play no role for the conductance. Channels with a significant width and amplitude occur around -2 eV, 0.5 eV and 3 eV. The corresponding admittance given in Fig. 2a) shows some similarities with the analytical model: below a resonance with , the conductance rises and the susceptance becomes more negative, as associated with a capacitive system. This can most clearly be seen near the first resonance at eV. Also at 2 eV and 3 eV, we can see further maxima in the real part of the admittance. A further comparison of numerical and analytical results is difficult, as the admittance of a device with multiple transmission channels will very likely differ from a simple superposition of single-level admittance traces. We are therefore currently exploring the extension of the Fu and Dudley result to the multi-level case.
4 Higher harmonics
Given a time-dependent bias of the form
one expects for non-linear devices a current response at multiples of the driving frequency . This can most easily be seen if one assumes that the time-dependent current follows the bias instantaneously and is well characterized by the current-voltage characteristic in the absence of harmonic driving. This allows one to write
Taking the expansion up to second order one obtains
These expressions explicitly show that a non-linear I-V will lead to higher harmonics in the current. Moreover, the time-averaged current will in general differ from the dc one, i.e., one observes current rectification. In recent experiments on nanoscale devices, the harmonic driving is realized by irradiation of the junction with frequencies in the microwave or optical range. Since the current is determined by the local voltage difference across the device, the ac bias entering Eq. 11 can not be obtained from the external field alone. Instead, the screening properties of the leads have to be taken into account, including especially also plasmonic excitations for higher frequencies. Assuming the validity of Eq. 12, this drawback might be turned into an advantage. Measuring both the photocurrent and , the local field enhancement may be estimated, which is difficult to access otherwise (39); (40); (22); (23).
Given this practical importance it is worthwhile to investigate under which circumstances Eqs. 11 and 12 are reasonable approximations. For Eq. 12 to hold, must be small and weakly non-linear around . For metallic nanoconstrictions the latter criterion is fulfilled for a relatively large bias range while for molecular junctions strong non-linearities arise close to molecular resonances. Around zero bias, the criterion is therefore more easily matched for fully saturated compounds compared to conjugated molecules, because of the larger gap in the former case.
The assumption on which the first line in Eq. 11 is based seems to be more severe. In a reactive circuit, current and voltage will in general be out of phase and it is not at all obvious that the dc I-V curve is relevant for the transport at high frequencies. In fact, the admittance calculations of the last section clearly showed that in molecular junctions already small frequency driving leads to sizable changes in the conductance. Moreover, Eq. 11 neglects quantum effects. In the Tien-Gordon picture for photo-assisted tunneling (24), the electrons may absorb multiple quanta of the radiation field in the leads before tunneling. As a consequence, the steady state transmission is probed at energies . Still, Eq. 12 might be a reasonable approximation in the limit of small , given that the transmission is weakly changing on the scale (41); (1).
In order to further investigate the limits of Eqs. 11 and 12, we now present atomistic calculations on the device depicted in Fig. 3, a 1,4-benzenediol molecule coupled to an Al nanowire. This junction was already characterized in earlier work (42); (35) and the admittance has also been computed in Ref. (38). Here we compute the higher harmonics in the current induced by an ac bias of the form given in Eq. 10. We generally take in order to ensure a smooth rise of the bias. This avoids high-frequency components in the current that are only related to initial transients and also increases the numerical stability of the time propagation. Both and can be independently adjusted in the experiment and were chosen to be on the same order in recent measurements (22); (23). Throughout the paper we work with an ac frequency of 200 THz, which correspond to a photon energy of 1.26 eV in the near-infrared.
We first compute the time-dependent current for low bias = 0.01 V with the method presented in Sec. 2. As shown in Fig. 4, the current leads
the voltage as expected for capacitive transport. Although the bias
polarity always remains positive, the current reverses sign. This
indicates that the conductance is larger than the dc
conductance. This matches with our earlier calculations on the
admittance of this device (38). Since the current may be evaluated
at the left and right device boundary, we also depict the sum of
both in Fig. 4. In the dc limit, and
coincide, while for a time-dependent bias the device region will in
general continuously charge and discharge. The displacement current we
find here is rather small
In order to examine the accuracy of the quasi-static picture of
Eq. 11 we proceed as follows. First, we compute the
curve using the conventional steady-state
Landauer-Büttiker formalism based on the DFTB Hamiltonian. As in Sec. 3, the same self energy is used in the static and
time-dependent calculations to ensure a common methodological
footing. The I-V curve is sampled in the interval [-2, 2] V with
an increment of 0.05 V and then fitted to a high-order polynomial. The
result is shown in Fig. 5a). Due to the
symmetric structure with equal lead-molecule distances
on both sides of the device, the I-V characteristic has
inversion symmetry around the origin. For low bias the current rises
roughly linearly, while around 2 V a steep increase is observed, which
arises due to the LUMO
Evaluating this function at the time-dependent bias Eq. 10, the quasi-static current in Fig. 5b) is obtained. Clearly, current and voltage are in phase and always of same polarity under this approximation in disagreement with the result given in Fig. 4. Further comparisons can be made by taking the Fourier transform of the current traces as given in Fig. 6. It can be seen that in the quasi-static approximation a larger number of harmonics exists, while the dynamic simulation shows a strong peak only at the fundamental frequency . As already evident from the polarity of the current in the time domain, we find . Since represents the dc component of the ac current, or in other words its time averaged value, we can also easily access the rectified current from the Fourier transform. Because of the small non-linearities in this low voltage regime (as also quantified in Fig. 6), we compute in the quasi-static limit a value of only 0.8 % for the ratio . Although the curvature at is positive (cf. Eq. 12), a negative value of -1.1 % is obtained for the dynamic simulations, showing that quantum effects can be important also for small dc bias.
Results for larger bias ( = 1 V) are given for the time domain in Fig. 7 and in the frequency domain in Fig. 8. While for the dynamic simulations the profile of I(t) is relatively smooth even at this larger bias, the quasi-static data has to follow the highly non-linear behaviour of around 2 V. Consequently, the corresponding Fourier transform features a larger number of higher harmonics. In contrast, only up to five harmonics are clearly discernible in the dynamic simulations which, however, have a larger amplitude. The rectified current reflects this with an increase of 60.8 % and 39.7 % for the dynamic and quasi-static simulations, respectively. Another difference between the two approaches is that the quasi-static falls off much stronger inbetween the harmonics. One possible reason for this could be that initial transients have not fully died out in the quantum calculations. This point could be further investigated by much longer simulations.
5 Conclusions and summary
In this article we reported on atomistic calculations of dynamical charge transport through single molecules. The employed time-domain approach based on density functional theory allows one to study arbitrary temporal profiles of the applied bias including in particular also ac driving. We find for molecular junctions that the classical quasi-static interpretation of high harmonic generation is of limited value already for small bias. If the ac frequency exceeds the microwave range, the conductance is significantly different from its dc value and the capacitance of the molecular device has to be taken into account, as discussed in Sec. 3. As a general trend, we observe that a full quantum treatment leads to smoother current traces and consequently also to a much smaller number of higher harmonics. For devices with a featureless transmission, the rectified current is proportional to dI/dV both in the classical quasi-static picture as well as in the quantum case (41); (1). Our results show that there is no such simple relation for molecular junctions. Dynamical quantum simulations should therefore play an important role in the interpretation of experiments on light driven ac transport in the future.
Acknowledgements.Financial support from the Deutsche Forschungsgemeinschaft (GRK 1570 and SPP 1243) is gratefully acknowledged.
- email: email@example.com
- email: firstname.lastname@example.org
- email: email@example.com
- journal: Journal of Computational Electronics
- For a deeper discussion of displacement currents see the introductory article by Oriols and Ferry in this special issue (43).
- LUMO = lowest unoccupied molecular orbital
- J.C. Cuevas, E. Scheer, Molecular Electronics: An Introduction to Theory and Experiment (World Scientific, 2010)
- M. Brandbyge, J.L. Mozos, P. Ordejon, J. Taylor, K. Stokbro, Phys. Rev. B 65(16), 165401 (2002)
- A.R. Rocha, V.M. García-Suárez, S. Bailey, C. Lambert, J. Ferrer, S. Sanvito, Phys. Rev. B 73(8), 085414 (2006)
- A. Pecchia, G. Penazzi, L. Salvucci, A. di Carlo, New J. Phys. 10, 065022 (2008)
- R.L. McCreery, H. Yan, A.J. Bergren, Phys. Chem. Chem. Phys. 15(4), 1065 (2013)
- D. Dulić, S. Van der Molen, T. Kudernac, H. Jonkman, J. De Jong, T. Bowden, J. Van Esch, B. Feringa, B. Van Wees, Phys. Rev. Lett. 91, 207402 (2003)
- S. van der Molen, J. Liao, T. Kudernac, J. Agustsson, L. Bernard, M. Calame, B. van Wees, B. Feringa, C. Schönenberger, Nano Letters 9, 76 (2008)
- S. Battacharyya, A. Kibel, G. Kodis, P.A. Liddell, M. Gervaldo, D. Gust, S. Lindsay, Nano Lett. 11(7), 2709 (2011)
- E.S. Tam, J.J. Parks, W.W. Shum, Y.W. Zhong, M.B. Santiago-Berrios, X. Zheng, W. Yang, G.K.L. Chan, H.D. Abruna, D.C. Ralph, Acs Nano 5(6), 5115 (2011)
- Y. Kim, T.J. Hellmuth, D. Sysoiev, F. Pauly, T. Pietsch, J. Wolf, A. Erbe, T. Huhn, U. Groth, U.E. Steiner, E. Scheer, Nano Lett. 12(7), 3736 (2012)
- X. Qiu, G. Nazin, W. Ho, Science 299(5606), 542 (2003)
- Y. Wakayama, K. Ogawa, T. Kubota, H. Suzuki, T. Kamikado, S. Mashiko, Applied Physics Letters 85, 329 (2004)
- Z.C. Dong, X.L. Guo, A. Trifonov, P. Dorozhkin, K. Miki, K. Kimura, S. Yokoyama, S. Mashiko, Phys. Rev. Lett. 92(8), 086801 (2004)
- Y. Fu, S. Dudley, Phys. Rev. Lett. 70(1), 65 (1993)
- Y. Yu, B. Wang, Y. Wei, J. Chem. Phys. 127(16), 169901 (2007)
- C.Y. Yam, Y. Mo, F. Wang, X. Li, G.H. Chen, X. Zheng, Y. Matsuda, J. Tahir-Kheli, A.G. William III, Nanotech. 19, 495203 (2008)
- T. Yamamoto, K. Sasaoka, S. Watanabe, Phys. Rev. B 82(20), 205404 (2010)
- T. Sasaoka, K.and Yamamoto, S. Watanabe, K. Shiraishi, Phys. Rev. B 84, 125403 (2011)
- D. Hirai, T. Yamamoto, S. Watanabe, Applied Physics Express 4(7), 075103 (2011)
- G.P. Kochanski, Phys. Rev. Lett. 62(19), 2285 (1989)
- S. Stranick, P. Weiss, J. Phys. Chem. 98(7), 1762 (1994)
- D.R. Ward, F. Hüser, F. Pauly, J.C. Cuevas, D. Natelson, Nature Nanotechnology 5(10), 732 (2010)
- R. Arielly, A. Ofarim, G. Noy, Y. Selzer, Nano Lett. 11(7), 2968 (2011)
- P.K. Tien, J.P. Gordon, Phys. Rev. 129, 647 (1963)
- J.K. Viljas, F. Pauly, J.C. Cuevas, Phys. Rev. B 76(3), 033403 (2007)
- J. Viljas, J. Cuevas, Phys. Rev. B 75(7), 075406 (2007)
- J.K. Viljas, F. Pauly, J.C. Cuevas, Phys. Rev. B 77(15), 155119 (2008)
- T. Brandes, Phys. Rev. B 56(3), 1213 (1997)
- G. Platero, R. Aguado, Physics Reports 395(1), 1 (2004)
- X. Zheng, F. Wang, C.Y. Yam, Y. Mo, G.H. Chen, Phys. Rev. B 75(19), 195127 (2007)
- X. Zheng, G.H. Chen, Y. Mo, S.K. Koo, H. Tian, C.Y. Yam, Y.J. Yan, J. Chem. Phys. 133, 114101 (2010)
- M.E. Casida, Recent Advances in Density Functional Methods, Part I (World Scientific, 1995), chap. Time-dependent Density Functional Response Theory for Molecules, pp. 155–192
- T. Frauenheim, G. Seifert, M. Elstner, T. Niehaus, C. Kohler, M. Amkreutz, M. Sternberg, Z. Hajnal, A. Di Carlo, S. Suhai, Journal Of Physics-Condensed Matter 14(11), 3015 (2002)
- T.A. Niehaus, J. Mol. Struct.: THEOCHEM 914, 38 (2009)
- Y. Wang, C.Y. Yam, T. Frauenheim, G. Chen, T. Niehaus, Chem. Phys. 391(1), 69 (2011)
- J. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77(18), 3865 (1996)
- B.N. Dutta, B. Dayal, physica status solidi (b) 3, 473 (1963)
- C. Oppenländer, B. Korff, T. Frauenheim, T.A. Niehaus, submitted to physica status solidi (b) arXiv:1304.4157 (2013). Submitted
- X. Tu, J. Lee, W. Ho, J. Chem. Phys. 124(2), 021105 (2006)
- D. Ward, G. Scott, Z. Keane, N. Halas, D. Natelson, J. Phys.: Condens. Matter 20, 374118 (2008)
- J.R. Tucker, M.J. Feldman, Rev. Mod. Phys. 57(4), 1055 (1985)
- C.Y. Yam, X. Zheng, G.H. Chen, Y. Wang, T. Frauenheim, T.A. Niehaus, Phys. Rev. B 83, 245448 (2011)
- X. Oriols, D. Ferry, J. Comp. Elec. (2013)