We present a time-dependent method based on the single-particle electron density matrix that allows the electronic and ionic degrees of freedom to be modelled within the Ehrenfest approximation in the presence of open boundaries. We describe a practical implementation using tight binding, and use it to investigate steady-state conduction through a single-atom device and to perform molecular dynamics. We find that in the Ehrenfest approximation an electric current allows both ionic heating and cooling to take place, depending on the bias

Bowler, D.R.

Fisher, A.J.

Horsfield, A.P.

English

UCL Discovery

LETTER TO THE EDITOROpen-boundary Ehrenfest molecular dynamics:towards a model of current induced heating innanowiresAndrew P. Hors¯eldE-mail: a.horsfield@ucl.ac.ukDepartment of Physics and Astronomy, University College London, GowerStreet, London WC1E 6BT, United KingdomD. R. BowlerE-mail: david.bowler@ucl.ac.ukDepartment of Physics and Astronomy, University College London, GowerStreet, London WC1E 6BT, United KingdomandLondon Centre for Nanotechnology, University College London, Gower St,London WC1E 6BTA. J. FisherE-mail: andrew.fisher@ucl.ac.ukDepartment of Physics and Astronomy, University College London, GowerStreet, London WC1E 6BT, United KingdomandLondon Centre for Nanotechnology, University College London, Gower St,London WC1E 6BTAbstract. We present a time-dependent method based on the single particleelectron density matrix that allows the electronic and ionic degrees of freedomto be modeled within the Ehrenfest approximation in the presence of openboundaries. We describe a practical implementation using tight binding, anduse it to investigate steady-state conduction through a single atom device and toperform molecular dynamics. We ¯nd that in the Ehrenfest approximation anelectric current allows both ionic heating and cooling to take place, depending onthe bias.Submitted to: J. Phys.: Condens. MatterLetter to the Editor 2One very important application of nanotechnology is to the development ofnanoscale electronic devices, possibly molecular[1, 2], that o®er higher performanceand lower power consumption than present technologies. As device sizes are scaleddown, heat production becomes an increasingly important design consideration.There are well developed modeling schemes now available to evaluate current-voltagecharacteristics for nanoscale devices[3, 4, 5, 6, 7, 8, 9, 10], but there are less satisfactoryschemes for evaluating the rate of heat generation[11].Within a Born-Oppenheimer framework, heating by an electric current occurs asa result of transitions between electronic states induced by atomic motion (includingzero point motion in the quantum case)[12, 11]. The prevailing atomistic modelsassume electrons permanently occupy well-de¯ned states characterised by two chemicalpotentials, and thus eliminate current induced heating, despite producing ionic forces.To progress we need an approach in which ions and electrons evolvesimultaneously in a consistent manner. There exist approaches of varyingcomplexity[13, 14, 15]. Here, we consider the simplest in which ions move along uniqueclassical trajectories determined by Hellmann-Feynman[16] forces (the Ehrenfestapproximation extended to open-boundary problems). We neglect all quantumcontributions to ionic motion. As we shall see, the Ehrenfest approximation is notadequate to model heating in all cases. However, in a paper in preparation we analyzethe reasons for the errors and propose a solution that is an extension of the methodpresented here.Here we describe a time-dependent formalism that is suitable for tight binding(TB) models implemented using density matrices. TB is chosen beause it is thesimplest quantum mechanical model of electrons that delivers quantitative results[17].Density matrices are used because they provide a very compact description of thestate of all the electrons[18, 19], and have proven very useful in the static descriptionof materials in the context of linear scaling methods[20]. Note that those parts of thedensity matrix treated explicitly must have a ¯nite range if they are to be used inpractical calculations. This restriction is elaborated on below.−−−−−+++++left right CRFigure 1. The circuit used to establish a current. The capacitor represents thenon-equilibrium source of charge, and the resistor the device through which wewish to drive the current.The physical model underpinning our method is represented in Fig. 1. A capacitorin series with a resistor form a complete circuit. For times t < 0, an external potential(equal to the chemical potential di®erence between the two plates of the capacitor)is applied to the left side of the circuit producing an excess of electrons on the leftand a de¯cit on the right. Most of the charge separation appears on the capacitor toLetter to the Editor 3minimise the energy.At time t = 0, the external potential is removed, and the electrons move so as toeliminate the imbalance, thereby producing a current °ow through the resistor and apotential drop across it. If ª0 is the many-body wave function for times t · 0 andthe many-body Hamiltonian for t ¸ 0 is ^ H, then the wave function for t ¸ 0 (ª(t))satis¯es ª(t) = exp( ^ Ht=i¹ h)ª0, provided @ ^ H=@t = 0. In the absence of dissipationthis leads to oscillatory solutions. However, if RC >> t >> ¹ h=W (W is the rangeof eigenvalues of ^ H whose corresponding vectors contribute to ª0), there will be aquasi-steady-state. It this time range in which we are interested. (This approach issimilar to that developed recently by Baer and Neuhauser[21].)In the quasi-steady-state regime most of the potential drop should occur acrossthe resistor. This allows us to focus on the resistor, and to treat the capacitor andmost of the wire as an external charge source and sink, which can be modeled by openboundary conditions. This is similar in spirit to the Landauer approach[3] in whichonly the transmission coe±cient for the device is evaluated.To simplify the time-dependent problem we reduce the full many-body formalismto an e®ective single particle one by working with time-dependent density functionaltheory (DFT)[22]. The key equations are:n(~ r;t) =XnfnjÃn(~ r;t)j2^ HksÃn(~ r;t) = i¹ h@@tÃn(~ r;t)^ Hks = ^ T + ^ VeI + ^ VHa[n] + ^ Vxc[n]MId2~ RIdt2 = ¡~ rIVII ¡Zd~ r n(~ r;t)~ rIVeI(~ r) (1)The last equation constitutes the Ehrenfest approximation. Here, n(~ r;t) is the chargedensity, fn is the orbital occupancy, Ãn(~ r;t) is a spin-dependent wavefunction of theKohn-Sham hamiltonian ^ Hks, ^ T is the kinetic energy operator, ^ VeI is the electron-ioninteraction, ^ VHa is the Hartree (electrostatic) interaction, ^ Vxc is the spin-dependentexchange and correlation potential, ^ VII is the ion-ion repulsion, MI is the mass of ion Iand ~ RI is its position. Note that charge self-consistency is automatically accomodatedby these equations. While it is perfectly possible to work directly with these equations(once suitable approximations for ^ Vxc have been made), we prefer to work with thesingle particle density matrix ½(~ r;~ r0;t), where½(~ r;~ r0;t) =XnÃn(~ r;t)fnÃ¤n(~ r0;t): (2)From Eqs (1) and (2) the equation of motion for the density matrices and the ionsare:i¹ h@^ ½@t= [ ^ Hks; ^ ½] (3)MId2~ RIdt2 = ¡~ rIVII ¡ Trf^ ½~ rI ^ Hksg (4)where we have moved to operator notation, for the density matrix.In this paper we use the TB approximation to DFT as this allows for analyticsimplicity and computational e±ciency[23, 24]. If we use orthogonal TB we canreplace operators in our previous equations with matrices. We thus continue to useLetter to the Editor 4the operator notation but with the new understanding that the operators will berepresented by TB matrices.We separate the system into the \resistor" (now refered to as the device) andthe remainder (the environment). We break Eq. (3) into components correspondingto device (designated by the subscript D), environment (designated by the subscriptE) and the coupling between the two. Further, we introduce a damping term for theenvironment:i¹ h@^ ½D@t= [ ^ HD; ^ ½D] + ( ^ HDE^ ½ED ¡ ^ ½DE ^ HED)i¹ h@^ ½DE@t= ^ HD^ ½DE ¡ ^ ½D ^ HDE + ^ HDE^ ½E ¡ ^ ½DE ^ HEi¹ h@^ ½E@t= [ ^ HE; ^ ½E] + ( ^ HED^ ½DE ¡ ^ ½ED ^ HDE)¡ 2i¹ h¡(^ ½E ¡ ^ ½ref): (5)The damping term corresponds to inelastic scattering events that take the system toa steady state con¯guration described by the reference density matrix ^ ½ref. In thesteady state, when we take the limit ¡ ! 0, we regain the key features of the quasi-static formalisms: the density matrix commutes with the Hamiltonian; the occupanciesof the single particle states are completely determined by the baths (represented by^ ½ref).The scattering rate ¡ is closely related to a self-energy[25]. It reduces the rangeof the density matrix in the environment, notably ^ ½ED which we must truncate tomake computations tractable. Truncating ^ ½ED without introducing damping resultsin the violation of charge conservation.We now ¯nd a closed form solution for the density matrix for the environment.We assume that ^ HE is independent of time (even in a self-consistent calculation thiscan be achieved by taking a su±ciently large device region) and de¯ne the driverterms ^ GE and ^ G(0)E by i¹ h ^ GE = ( ^ HED^ ½DE ¡ ^ ½ED ^ HDE) and 0 = [ ^ HE; ^ ½E(0)]+i¹ h ^ G(0)E ¡2i¹ h¡(^ ½E(0) ¡ ^ ½ref). By making use of the interaction picture the solution for Eq. (5)for the environment is found to be^ ½E(t) = ^ ½E(0) +Z t0dx ^ O(x)³^ GE(t ¡ x) ¡ ^ G(0)E´^ Oy(x) (6)where ^ O(t) = e¡¡te^ HEt=i¹ h. The time evolution matrices, ^ O(t), are moststraightforwardly evaluated from Green's functions:^ O(t) = e¡¡tZdE eEt=i¹ h±(E^ IE ¡ ^ HE) (7)where±(E^ IE ¡ ^ HE) = ^ d(E) = ¡1¼lim´!0+ Im ^ g(E + i´): (8)provided ^ HE is real. Here ^ g(Z) is the one particle Green's function for the environmentwhich can be evaluated using standard methods[26, 27].We cannot use the closed form solution for the part of the density matrixbelonging to the device and its coupling to the environment as the relevant partsof the Hamiltonian will vary with time due to the atomic motion we wish to study.Instead we treat the time evolution explicitly. We can do this because this subsystemLetter to the Editor 5is small, provided the density matrix ^ ½DE is short ranged. We evolve the relevantparts of the density matrix using the approximation f(t+¢t) ¼ f(t¡¢t)+2¢tf0(t).Note that the above formalism is similar in spirit to a recent wave functionbased method[28]. However, our use of the density matrix allows for a more naturalinterface between the device and the environment as we do not have to treat individualwavefunctions but only an integrated quantity.To see how the method behaves, consider a very simple model system. It consistsof two semi-in¯nite leads attached to a device. The lead on the left is at a di®erentpotential from that on the right. Each lead is represented by a linear chain of atomswith one orbital per atom. Thus there are two parameters that characterise theHamiltonian for each lead: an on-site energy (a) and a hopping integral between thenearest neighbour sites (b). We assume b is the same on the left and right, while thedi®erence in onsite energies aL ¡ aR corresponds to the potential di®erence betweenthe two sides. We take both the initial and reference density matrices (^ ½E(0) and^ ½ref) to be that for the in¯nite wire with the device in its ground state in the absenceof a bias.If we label the atoms in one lead 0;1;2:::, where 0 is next to the device, then wecan de¯ne the time evolution matrix (see Eqs. (7) and (8)):Onm(t) = e¡¡tZ a+2jbja¡2jbjdE eEt=i¹ hdnm(E): (9)This integral is evaluated numerically.0 2 4 6Time (fs)00.51O00Figure 2. The solid line is the variation with time of the matrix element ofthe time evolution operator corresponding to the ¯rst atom in the environment(O00(t)). There is no damping (¡ = 0). The decay corresponds to the propagationof a wavepacket down the wire. The dashed line is the exponential damping factorwith ¡ = 1:0fs¡1.The non-locality in time of Eq. (6) adds considerably to the cost of performing acalculation. However, from Fig. 2 we see that the evolution operator decays rapidlywith time. We thus approximate this by a function that goes to zero after a cut-o®time. This is consistent with keeping a ¯nite ¡ in Eq. (7). If the evolution operator istruncated in time it is truncated in space as well.The simplest device consists of one atom. If we give this atom a high on-siteenergy, it behaves as a barrier to current °ow. The one-dimensional potential pro¯lefor this system is given in Fig. 3. As a check on the method described here, wecomputed the conductivity of this system using the Landauer method[3]. With a biasof 0.1V, and a hopping integral and barrier height both of 1eV, it gives a currentof about 6.2¹A. It should be noted, however, that there is no guarantee that theLetter to the Editor 6DeviceaaaRLDLeft leadRight leadEnergyFigure 3. The energy pro¯le for the model system. The energy axis on the leftshows the positions of the onsite energies in the left lead (aL), the right lead (aR)and the device (aD).time-dependent formalism will exactly reproduce the Landauer result. The orbitaloccupancies in the Landauer theory are de¯ned by two chemical potentials (¯xedby the bias and charge neutrality). In the present theory they depend on the bias,reference density matrix, and ¡ if it is large.The time dependence of the current is shown in Fig. 4(a). We see that our schemeleads to stable steady currents. We have tested convergence of our results with devicesize and the range of the density matrix from the device into the leads. For the devicewe allow only one atom to have the increased onsite energy with the remaining atomsbeing lead atoms. Results are summarised in Fig. 5. The current is generally belowthe Landauer value and is insensitive to the range of the density matrix from thedevice into the leads. As the reference density matrix is real, the imaginary part of-20246810Current (µA)0 20 40 60 80 100Time (fs)-202468Current (µA)(a)(b)Figure 4. The variation of current into the device as a function of time. Notethat it reaches a stable steady state. For the ¯rst 30 fs no bias was applied. Thebias was then turned on over a period of 10 fs. The top panel (a) shows thesystem with no self-consistency applied, the bottom (b) with self consistency, asdescribed in the text.Letter to the Editor 75 10 15 20Number of atoms in device55.566.5Current (µΑ)Landauer8 atom lead16 atom leadFigure 5. The steady currents calculated as a function of the system size. Noteonly one atom in the middle has the increased onsite energy of 1eV. The bias is0.1V, the hopping integral is 1eV and ¹ h¡ = 0:0658 eV.the density matrix is suppressed by the damping term, thereby reducing the currentin the leads. Damping dominates the electron dynamics when ¹ h¡=W ¸ 1, where Wis the band width. In the present case ¹ h¡=W = 0:016.We have repeated a conduction calculation with charge self-consistency included.We use a simple point monopole scheme: if the net charge on site i is Qi, then itcontributes a potential UiQi (Ui=7eV) to its own site and Qi=qR2ij + 1=U2i (whereRij and Ui are in atomic units) to neighbouring sites. To match the potential at theends of the device to that in the leads we include a linear external ¯eld in the deviceregion of V (x) = a+bx. The results are shown in Fig. 4(b). It can be seen that, asidefrom a sunstantial increase in the noise in the system (which decays away rapidly),there is little change in the overall behaviour.We now come to the calculation for which the method is designed. Once a steady-state current has been achieved, we assign a velocity to the device atom, correspondingto a given temperature, and then perform molecular dynamics on it in one dimension.The position of the ion is allowed to evolve according to Eq. (4). Note the environmentatoms do not move, so no heat can be transported away.0.000.050.100.15K.E. (eV) 0 100 200 300 400Time (fs)0.000.050.10K.E. (eV)a)b)Figure 6. The variation of the ionic kinetic energy of the mobile atom with timefor a bias of a) 0.1V and b) 1.0V. The device contains 3 atoms, and the leads 16atoms. The initial temperature is about 600K.Letter to the Editor 8To monitor ionic heating we follow the evolution of the kinetic energy with time(see Fig. 6). For a small bias (0.1V) the ionic kinetic energy decays with time (cooling),while for a large bias (1.0V) the kinetic energy increases (heating). The energy (¹ h!)associated with the vibrational frequency of the ion is 0.055 eV, and equals the Born-Oppenheimer surface separation for allowed transitions. Hence the change in biasincreases approximately tenfold the number of possible heating transitions, producingthe observed changed behaviour.This method is thus able to calculate some non-adiabatic e®ects of current °ow.However, the amount of heating that we observe for a bias of 1V is much less than wewould expect from quantum perturbation theory[12]. The explanation for this, anda correction to remedy it, are the subject of ongoing work and will be presented in apaper currently in preparation.AcknowledgmentsThe ¯nancial support of the IRC on nanotechnology for APH and of the Royal Societyfor DRB, and conversations with Tchavdar Todorov, Adrian Sutton and MarshallStoneham are gratefully acknowledged.[1] Reed M A, Zhou C, Muller C J, Burgin T P, and Tour J M 1997 Science 278(5336) 252[2] Di Ventra M, Pantelides S T, and Lang N D 2000 Appl. Phys. Lett. 76(23) 3448[3] Landauer R 1989 J. Phys.: Condens. Matter 1 8099[4] Lang N 1995 Phys. Rev. B 52(7) 5335[5] Ness H and Fisher A J 1997 Phys. Rev. B 56(19) 12469[6] Sanvito S, Lambert C J, Je®erson J H, and Bratkovsky A M 1999 Phys. Rev. B 59(18) 11936[7] Taylor J, Guo H, and Wang J 2001 Phys. Rev. B 63(24) 245407[8] Brandbyge M, Mozos J L, Ordejon P, Taylor J, and Stokbro K 2002 Phys. Rev. B 65(16) 165401[9] Todorov T N 2002 J. Phys.: Condens. Matter 14 3049[10] Di Ventra M and Lang N D 2002 Phys. Rev. B 65(4) 45402[11] Montgomery M J, Todorov T N, and Sutton A P 2002 J. Phys.: Condens. Matter 14 5377[12] Todorov T N 1998 Phil. Mag. B 77(4) 965[13] Kohen D, Stillinger F H, and Tully J C 1998 J. Chem. Phys. 109(12) 4713[14] Billing G D 1999 J. Chem. Phys. 110(12) 5526[15] Kapral R and Ciccotti G 1999 J. Chem. Phys. 110(18) 8919[16] Di Ventra M and Pantelides S T 2000 Phys. Rev. B 61(23) 16207[17] Goringe C M, Bowler D R, and Hern¶ andez E 1997 Rep. Prog. Phys. 60(12) 1447[18] McWeeny R 1960 Rev. Mod. Phys. 32 335[19] He L and Vanderbilt D 2001 Phys. Rev. Lett. 86(23) 5341[20] Goedecker S 1999 Rev. Mod. Phys. 71(4) 1085[21] Baer R and Neuhauser D 2003 Chem. Phys. Lett. 374 459[22] Calvayrac F, Reinhard P G, Suraud E, and Ullrich C A 2000 Physics Reports 337(6) 493[23] Hors¯eld A P and Bratkovsky A M 2002 J. Phys.: Condens. Matter 12 R1[24] Todorov T N 2001 J. Phys.: Condens. Matter 13 10125[25] Todorov T N 2003 Private communication[26] Williams A R, Feibelman P J, and Lang N D 1982 Phys. Rev. B 26(10) 5433[27] L¶ opez Sancho M P, L¶ opez Sancho J M, and Rubio J 1984 J. Phys. F: Met. Phys. 14 1205[28] Baer R and Neuhauser D 2003 Int. J. Quantum Chem. 91 524

Open-boundary Ehrenfest molecular dynamics: towards a model of current unduced heating in nanowires

Abstract

Similar works

Full text

Available Versions

UCL Discovery