- Email: [email protected]
Self-consistent conductance calculations on molecular calipers using a transfer matrix method
Superlattices and Microstructures 34 (2003) 429–432 www.elsevier.com/locate/superlattices
Self-consistent conductance calculations on molecular calipers using a transfer matrix method G. Speyer∗, R. Akis, D.K. Ferry Department of Electrical Engineering and Center for Solid State Electronics Research, Arizona State University, Tempe, AZ 85287-5706, USA Available online 15 June 2004
Abstract A rapid method of conductance calculation using an iterative transfer matrix algorithm has been developed. Its use is demonstrated with a simple metal–molecule–metal system. By using a DFT Hamiltonian and charge self-consistency, we avoid any dependence on fitting. Here we solve a selfconsistent potential, which obviates the need to parametrize the voltage. Moreover, in examining an experimental set-up developed at Arizona State University, the molecular conductance across a variety of gap lengths can be calculated and compared to experiment, and the implementation of a proposed nano-caliper device can be evaluated. Conduction across the molecule occurs in multiple channels; gold states couple with varying strengths to the HOMO of the molecule. We will report the effects of strain across the molecule, and of distortion of the molecule, on the conductive nature of the coupling. © 2004 Elsevier Ltd. All rights reserved. Keywords: Molecular electronics; Transfer matrix method
Conductance calculations to this point have either followed direct solution of scattering states [1] or employed non-equilibrium Green’s functions (NEGF) [2]. The former involves explicit solution for all the eigenvectors of a large rank matrix to calculate the wavefunctions for density. The NEGF approach calculates the transmission matrix at a specific electron energy and bias voltage by taking the integrated trace of the product of the contact-dependent gammas (the imaginary part of the contact self-energy) and the molecule-specific Green’s function matrices. The surface Green’s functions used in NEGF must be calculated using a recursive minimization technique at each energy point used in the integration. In order to calculate the density, the molecular Green’s function and the ∗ Corresponding author. Tel.: +1-480-965-3452.
E-mail address: [email protected] (G. Speyer). 0749-6036/$ - see front matter © 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.spmi.2004.03.063
430
G. Speyer et al. / Superlattices and Microstructures 34 (2003) 429–432
Fig. 1. Discretization of the Hamiltonian onto the lattice. Molecular sites are consolidated onto a single slice in order to maximize the number of nearest neighbors.
correlation matrix must be integrated to yield the equilibrium and non-equilibrium matrices respectively. Due to the fine structure of the density of states along the real axis, these integrals must be carefully calculated. Different approaches employ analytic continuation which, in turn, implies the surface Green’s functions must also be calculated along the complex contour. Even if the surface Green’s function calculation, which involves several matrix inversions among other time-intensive operations, is calculated a priori, the amount of calculation becomes cumbersome. Our approach exploits the computational efficiency and numerical stability of a transfer matrix variant. By discretizing the system into sites on a lattice (see Fig. 1), the transmission is solved through iterative matrix calculation and translation at each slice, much like in a recursive Green’s function technique. Starting with Hamiltonian matrix elements and interatomic distances for the system, the molecular and contact atomic sites are mapped onto a discrete lattice whereby an incident electron flux is propagated slice by slice. The transmission coefficients yielded at the last slice enter the Landauer–Buttiker formula to give conductance [3]. The divergence factors which contribute to the instability of traditional transfer matrix methods are shrewdly avoided using matrix manipulations which cancel these factors. The advantage over Green’s function methods is that the transmission matrix and wavefunction for all points can be derived simultaneously. This saves considerable time in the self-consistency scheme in which the density (i.e., the normalized square of the wavefunction) must be repetitively calculated. Self-consistency is achieved using a direct integration of a pre-calculated Green’s function, G(r, r0 ), multiplied by the charge terms, ρ(r ), at each site with Dirichlet boundary conditions, f (r ), imposed: ∂ G(r, r0 ) u(r0 ) = dS(r ). (1) G(r, r0 )ρ(r ) dV (r ) + f (r ) ∂n V S In order to solve the Dirichlet problem, the Green’s function must have zero value at the boundaries.1 This is accomplished by creating image sources of opposite polarity.
1 NB: These Green’s functions are simply used to handle the electrostatic interactions. These are not the Green’s functions we seek to avoid using, namely those which are used to calculate transmission through the coupled system.
G. Speyer et al. / Superlattices and Microstructures 34 (2003) 429–432
431
Fig. 2. Xylyldithiol in a variable gap between two gold contacts.
However, due to there being two boundaries, the images of the images must be created as well. One-dimensional Ewald sums over these periodic arrays of charges are therefore implemented for every charge inside the system [4]. The first-principles program FIREBALL 2000, a local atomic orbital density functional theory (DFT) based method in the pseudopotential local density approximation (LDA), was used to calculate the Hamiltonian employed by the transfer matrix code [5]. Some previous investigations into this molecular system employed either an empirical tightbinding or an extended H¨uckel model and parametrized the voltage division factor, phonon coupling strength, and dielectric constant [6]. Fig. 2 shows the unit cell. In addition, FIREBALL 2000 was also used to calculate the Hellman–Feynman forces used to drift the xylyldithiol atoms upon stretching. In order to preserve the periodicity of the unit cell, the gold atoms were left fixed in these simulations, although the dynamics of the gold atoms are believed to be important in the stretching. Molecules were initially attached in the hollow-site configuration. The experimental system being modeled was developed by Tao and co-workers at Arizona State University and involves the formation of thousands of molecular junctions for measurement by repeatedly lowering a gold-plated AFM tip into and out of a gold substrate where a solution of xylyldithiol molecules has been deposited [7]. Because each molecule is terminated by thiol groups, the molecular junction readily forms between tip and substrate (Fig. 2), and indeed traces indicating conductance through discrete numbers of molecules were measured. Our calculations agree within an order of magnitude with experimental calculations of the molecular conductance and indicate an interesting trend in the conductance upon ˚ Orbital plots (Fig. 3(b) and (c)) stretching (Fig. 3(a)), with an apparent resonance for ∼2 A. help explain this phenomenon. As the molecule is stretched, orbitals near the Fermi level change in degree of localization. At the resonance, there is a conductance enhancement due to the planarization of molecule leading to enhanced coupling between gold states and molecular states. This is evident in the LUMO-like orbital shown in Fig. 3(c). We also note
432
G. Speyer et al. / Superlattices and Microstructures 34 (2003) 429–432
Fig. 3. (a) The conductance trend as the molecule is stretched. The spike reflects a similar phenomenon observed ˚ stretching are shown in (b) and (c). Planarization experimentally. The LUMO-like wavefunctions for 0.6 and 2 A of the molecule enhances the orbital delocalization.
the effect of charge transfer at the metal–molecule interface with applied bias, agreeing well with other theoretical observations [8]. Acknowledgements This work was supported by the Office of Naval Research. References [1] M. Di Ventra, N.D. Lang, Phys. Rev. B 65 (2001) 454021; E. Emberly, G. Kirczenow, Phys. Rev. B 58 (1998) 10911. [2] S. Datta, Electronic Transport in Mesoscopic Systems, Cambridge University Press, Cambridge, 1997; V. Mujica, A. Roitberg, M. Ratner, J. Chem. Phys. 112 (15) (2000) 6834; The similarity of these Green’s Functions techniques is explained in L. Hall, J.R. Reimers, N.S. Hush, K. Silverbrook, J. Chem. Phys. 112 (3) (2000) 1510. [3] T. Usuki, M. Saito, M. Takatsu, R.A. Kiehl, N. Yokoyama, Phys. Rev. B 52 (1995) 8244. [4] D.J. Langridge, J.F. Hart, S. Crampin, Comput. Phys. Comm. 134 (2001) 78. [5] O.F. Sankey, D.J. Niklewski, Phys. Rev. B 40 (1989) 3979. [6] W. Tian, S. Datta, S. Hong, R. Reifenberger, J. Henderson, C. Kubiak, J. Chem. Phys. 109 (1998) 2874; P. Kornilovitch, A. Bratkovsky, Phys. Rev. B 64 (2001) 195413; A. Tikhonov, R.D. Coalson, Y. Dahnovsky, J. Chem. Phys. 117 (2002) 567. [7] B. Xu, N.J. Tao, Science 301 (2003) 1221. [8] Y. Xue, M.A. Ratner, Phys. Rev. B 68 (2003) 115406.
Self-consistent conductance calculations on molecular calipers using a transfer matrix method G. Speyer∗, R. Akis, D.K. Ferry Department of Electrical Engineering and Center for Solid State Electronics Research, Arizona State University, Tempe, AZ 85287-5706, USA Available online 15 June 2004
Abstract A rapid method of conductance calculation using an iterative transfer matrix algorithm has been developed. Its use is demonstrated with a simple metal–molecule–metal system. By using a DFT Hamiltonian and charge self-consistency, we avoid any dependence on fitting. Here we solve a selfconsistent potential, which obviates the need to parametrize the voltage. Moreover, in examining an experimental set-up developed at Arizona State University, the molecular conductance across a variety of gap lengths can be calculated and compared to experiment, and the implementation of a proposed nano-caliper device can be evaluated. Conduction across the molecule occurs in multiple channels; gold states couple with varying strengths to the HOMO of the molecule. We will report the effects of strain across the molecule, and of distortion of the molecule, on the conductive nature of the coupling. © 2004 Elsevier Ltd. All rights reserved. Keywords: Molecular electronics; Transfer matrix method
Conductance calculations to this point have either followed direct solution of scattering states [1] or employed non-equilibrium Green’s functions (NEGF) [2]. The former involves explicit solution for all the eigenvectors of a large rank matrix to calculate the wavefunctions for density. The NEGF approach calculates the transmission matrix at a specific electron energy and bias voltage by taking the integrated trace of the product of the contact-dependent gammas (the imaginary part of the contact self-energy) and the molecule-specific Green’s function matrices. The surface Green’s functions used in NEGF must be calculated using a recursive minimization technique at each energy point used in the integration. In order to calculate the density, the molecular Green’s function and the ∗ Corresponding author. Tel.: +1-480-965-3452.
E-mail address: [email protected] (G. Speyer). 0749-6036/$ - see front matter © 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.spmi.2004.03.063
430
G. Speyer et al. / Superlattices and Microstructures 34 (2003) 429–432
Fig. 1. Discretization of the Hamiltonian onto the lattice. Molecular sites are consolidated onto a single slice in order to maximize the number of nearest neighbors.
correlation matrix must be integrated to yield the equilibrium and non-equilibrium matrices respectively. Due to the fine structure of the density of states along the real axis, these integrals must be carefully calculated. Different approaches employ analytic continuation which, in turn, implies the surface Green’s functions must also be calculated along the complex contour. Even if the surface Green’s function calculation, which involves several matrix inversions among other time-intensive operations, is calculated a priori, the amount of calculation becomes cumbersome. Our approach exploits the computational efficiency and numerical stability of a transfer matrix variant. By discretizing the system into sites on a lattice (see Fig. 1), the transmission is solved through iterative matrix calculation and translation at each slice, much like in a recursive Green’s function technique. Starting with Hamiltonian matrix elements and interatomic distances for the system, the molecular and contact atomic sites are mapped onto a discrete lattice whereby an incident electron flux is propagated slice by slice. The transmission coefficients yielded at the last slice enter the Landauer–Buttiker formula to give conductance [3]. The divergence factors which contribute to the instability of traditional transfer matrix methods are shrewdly avoided using matrix manipulations which cancel these factors. The advantage over Green’s function methods is that the transmission matrix and wavefunction for all points can be derived simultaneously. This saves considerable time in the self-consistency scheme in which the density (i.e., the normalized square of the wavefunction) must be repetitively calculated. Self-consistency is achieved using a direct integration of a pre-calculated Green’s function, G(r, r0 ), multiplied by the charge terms, ρ(r ), at each site with Dirichlet boundary conditions, f (r ), imposed: ∂ G(r, r0 ) u(r0 ) = dS(r ). (1) G(r, r0 )ρ(r ) dV (r ) + f (r ) ∂n V S In order to solve the Dirichlet problem, the Green’s function must have zero value at the boundaries.1 This is accomplished by creating image sources of opposite polarity.
1 NB: These Green’s functions are simply used to handle the electrostatic interactions. These are not the Green’s functions we seek to avoid using, namely those which are used to calculate transmission through the coupled system.
G. Speyer et al. / Superlattices and Microstructures 34 (2003) 429–432
431
Fig. 2. Xylyldithiol in a variable gap between two gold contacts.
However, due to there being two boundaries, the images of the images must be created as well. One-dimensional Ewald sums over these periodic arrays of charges are therefore implemented for every charge inside the system [4]. The first-principles program FIREBALL 2000, a local atomic orbital density functional theory (DFT) based method in the pseudopotential local density approximation (LDA), was used to calculate the Hamiltonian employed by the transfer matrix code [5]. Some previous investigations into this molecular system employed either an empirical tightbinding or an extended H¨uckel model and parametrized the voltage division factor, phonon coupling strength, and dielectric constant [6]. Fig. 2 shows the unit cell. In addition, FIREBALL 2000 was also used to calculate the Hellman–Feynman forces used to drift the xylyldithiol atoms upon stretching. In order to preserve the periodicity of the unit cell, the gold atoms were left fixed in these simulations, although the dynamics of the gold atoms are believed to be important in the stretching. Molecules were initially attached in the hollow-site configuration. The experimental system being modeled was developed by Tao and co-workers at Arizona State University and involves the formation of thousands of molecular junctions for measurement by repeatedly lowering a gold-plated AFM tip into and out of a gold substrate where a solution of xylyldithiol molecules has been deposited [7]. Because each molecule is terminated by thiol groups, the molecular junction readily forms between tip and substrate (Fig. 2), and indeed traces indicating conductance through discrete numbers of molecules were measured. Our calculations agree within an order of magnitude with experimental calculations of the molecular conductance and indicate an interesting trend in the conductance upon ˚ Orbital plots (Fig. 3(b) and (c)) stretching (Fig. 3(a)), with an apparent resonance for ∼2 A. help explain this phenomenon. As the molecule is stretched, orbitals near the Fermi level change in degree of localization. At the resonance, there is a conductance enhancement due to the planarization of molecule leading to enhanced coupling between gold states and molecular states. This is evident in the LUMO-like orbital shown in Fig. 3(c). We also note
432
G. Speyer et al. / Superlattices and Microstructures 34 (2003) 429–432
Fig. 3. (a) The conductance trend as the molecule is stretched. The spike reflects a similar phenomenon observed ˚ stretching are shown in (b) and (c). Planarization experimentally. The LUMO-like wavefunctions for 0.6 and 2 A of the molecule enhances the orbital delocalization.
the effect of charge transfer at the metal–molecule interface with applied bias, agreeing well with other theoretical observations [8]. Acknowledgements This work was supported by the Office of Naval Research. References [1] M. Di Ventra, N.D. Lang, Phys. Rev. B 65 (2001) 454021; E. Emberly, G. Kirczenow, Phys. Rev. B 58 (1998) 10911. [2] S. Datta, Electronic Transport in Mesoscopic Systems, Cambridge University Press, Cambridge, 1997; V. Mujica, A. Roitberg, M. Ratner, J. Chem. Phys. 112 (15) (2000) 6834; The similarity of these Green’s Functions techniques is explained in L. Hall, J.R. Reimers, N.S. Hush, K. Silverbrook, J. Chem. Phys. 112 (3) (2000) 1510. [3] T. Usuki, M. Saito, M. Takatsu, R.A. Kiehl, N. Yokoyama, Phys. Rev. B 52 (1995) 8244. [4] D.J. Langridge, J.F. Hart, S. Crampin, Comput. Phys. Comm. 134 (2001) 78. [5] O.F. Sankey, D.J. Niklewski, Phys. Rev. B 40 (1989) 3979. [6] W. Tian, S. Datta, S. Hong, R. Reifenberger, J. Henderson, C. Kubiak, J. Chem. Phys. 109 (1998) 2874; P. Kornilovitch, A. Bratkovsky, Phys. Rev. B 64 (2001) 195413; A. Tikhonov, R.D. Coalson, Y. Dahnovsky, J. Chem. Phys. 117 (2002) 567. [7] B. Xu, N.J. Tao, Science 301 (2003) 1221. [8] Y. Xue, M.A. Ratner, Phys. Rev. B 68 (2003) 115406.