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

Andrew P Horsfield¹, D R Bowler^1,2 and A J Fisher^1,2

¹Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK
²London Centre for Nanotechnology, University College London, Gower Street, London WC1E 6BT, UK

Email: a.horsfield@ucl.ac.uk, david.bowler@ucl.ac.uk and andrew.fisher@ucl.ac.uk

Received 22 January 2004
Published 6 February 2004

PACS number: 7363N

One very important application of nanotechnology is to the development of nanoscale electronic devices, possibly molecular [1, 2], that offer higher performance and lower power consumption than present technologies. As device sizes are scaled down, heat production becomes an increasingly important design consideration. There are well-developed modelling schemes now available to evaluate current-voltage characteristics for nanoscale devices [3-10], but there are fewer satisfactory schemes for evaluating the rate of heat generation [11].

Within a Born-Oppenheimer framework, heating by an electric current occurs as a result of transitions between electronic states induced by atomic motion (including zero-point motion in the quantum case) [12, 11]. The prevailing atomistic models assume that electrons permanently occupy well-defined states characterized by two chemical potentials, and thus eliminate current induced heating, despite producing ionic forces.

To progress we need an approach in which ions and electrons evolve simultaneously in a consistent manner. There exist approaches of varying complexity [13-15]. Here, we consider the simplest, in which ions move along unique classical trajectories determined by Hellmann-Feynman [16] forces (the Ehrenfest approximation extended to open-boundary problems). We neglect all quantum contributions to ionic motion. As we shall see, the Ehrenfest approximation is not adequate to model heating in all cases. However, in a paper in preparation we analyse the reasons for the errors and propose a solution that is an extension of the method presented here.

Here we describe a time-dependent formalism that is suitable for tight binding (TB) models implemented using density matrices. TB is chosen because it is the simplest quantum mechanical model of electrons that delivers quantitative results [17]. Density matrices are used because they provide a very compact description of the state of all the electrons [18, 19], and have proven very useful in the static description of materials in the context of linear scaling methods [20]. Note that those parts of the density matrix treated explicitly must have a finite range if they are to be used in practical calculations. This restriction is elaborated on below.

The physical model underpinning our method is represented in figure 1. A capacitor in series with a resistor forms a complete circuit. For times t<0, an external potential (equal to the chemical potential difference between the two plates of the capacitor) is applied to the left side of the circuit producing an excess of electrons on the left and a deficit on the right. Most of the charge separation appears on the capacitor to minimize the energy.

At time t  =  0, the external potential is removed, and the electrons move so as to eliminate the imbalance, thereby producing a current flow through the resistor and a potential drop across it. If Ψ₀ is the many-body wavefunction for times and the many-body Hamiltonian for is , then the wavefunction for (Ψ(t)) satisfies , provided that . In the absence of dissipation this leads to oscillatory solutions. However, if (W is the range of eigenvalues of whose corresponding vectors contribute to Ψ₀), there will be a quasi-steady-state. It is this time range in which we are interested. (This approach is similar to that developed recently by Baer and Neuhauser [21].)

In the quasi-steady-state regime most of the potential drop should occur across the resistor. This allows us to focus on the resistor, and to treat the capacitor and most of the wire as an external charge source and sink, which can be modelled by open boundary conditions. This is similar in spirit to the Landauer approach [3] in which only the transmission coefficient for the device is evaluated.

To simplify the time-dependent problem we reduce the full many-body formalism to an effective single-particle one by working with time-dependent density functional theory (DFT) [22]. The key equations are

The last equation constitutes the Ehrenfest approximation. Here, is the charge density, f_n is the orbital occupancy, is a spin-dependent wavefunction of the Kohn-Sham Hamiltonian , is the kinetic energy operator, is the electron-ion interaction, is the Hartree (electrostatic) interaction, is the spin-dependent exchange and correlation potential, is the ion-ion repulsion, M_I is the mass of ion I and is its position. Note that charge self-consistency is automatically accommodated by these equations. While it is perfectly possible to work directly with these equations (once suitable approximations for have been made), we prefer to work with the single-particle density matrix , where

From equations (1) and (2) the equations of motion for the density matrices and the ions are

where we have moved to operator notation for the density matrix.

In this letter we use the TB approximation to DFT as this allows for analytic simplicity and computational efficiency [23, 24]. If we use orthogonal TB we can replace operators in our previous equations with matrices. We thus continue to use the operator notation but with the new understanding that the operators will be represented by TB matrices.

We separate the system into the `resistor' (now referred to as the device) and the remainder (the environment). We break equation (3) into components corresponding to the device (designated by the subscript D), the environment (designated by the subscript E) and the coupling between the two. Further, we introduce a damping term for the environment:

The damping term corresponds to inelastic scattering events that take the system to a steady-state configuration described by the reference density matrix . In the steady state, when we take the limit , we regain the key features of the quasi-static formalisms: the density matrix commutes with the Hamiltonian; the occupancies of the single-particle states are completely determined by the baths (represented by ).

The scattering rate Γ is closely related to a self-energy [25]. It reduces the range of the density matrix in the environment, notably , which we must truncate to make computations tractable. Truncating  without introducing damping results in the violation of charge conservation.

We now find a closed-form solution for the density matrix for the environment. We assume that is independent of time (even in a self-consistent calculation this can be achieved by taking a sufficiently large device region) and define the driver terms and by and . By making use of the interaction picture, the solution for equation (5) for the environment is found to be

where . The time evolution matrices, , are most straightforwardly evaluated from Green functions:

where

provided that is real. Here is the one-particle Green function for the environment which can be evaluated using standard methods [26, 27].

We cannot use the closed-form solution for the part of the density matrix belonging to the device and its coupling to the environment, as the relevant parts of the Hamiltonian will vary with time due to the atomic motion that we wish to study. Instead we treat the time evolution explicitly. We can do this because this subsystem is small, provided that the density matrix is short ranged. We evolve the relevant parts of the density matrix using the approximation .

Note that the above formalism is similar in spirit to a recent wavefunction based method [28]. However, our use of the density matrix allows for a more natural interface between the device and the environment as we do not have to treat individual wavefunctions but only an integrated quantity.

To see how the method behaves, consider a very simple model system. It consists of two semi-infinite leads attached to a device. The lead on the left is at a different potential to that on the right. Each lead is represented by a linear chain of atoms with one orbital per atom. Thus there are two parameters that characterize the Hamiltonian for each lead: an on-site energy (a) and a hopping integral between the nearest neighbour sites (b). We assume that b is the same on the left and right, while the difference in on-site energies a_L -  a_R corresponds to the potential difference between the two sides. We take both the initial and reference density matrices ( and ) to be that for the infinite wire with the device in its ground state in the absence of a bias.

If we label the atoms in one lead 0,1,2..., where 0 is next to the device, then we can define the time evolution matrix (see equations (7) and (8)):

This integral is evaluated numerically.

The non-locality in time of equation (6) adds considerably to the cost of performing a calculation. However, from figure 2 we see that the evolution operator decays rapidly with time. We thus approximate this by a function that goes to zero after a cut-off time. This is consistent with keeping a finite Γ in equation (7). If the evolution operator is truncated in time it is truncated in space as well.

The simplest device consists of one atom. If we give this atom a high on-site energy, it behaves as a barrier to current flow. The one-dimensional potential profile for this system is given in figure 3. As a check on the method described here, we computed the conductivity of this system using the Landauer method [3]. With a bias of 0.1 V, and a hopping integral and barrier height both of 1 eV, it gives a current of about 6.2 µA. It should be noted, however, that there is no guarantee that the time-dependent formalism will exactly reproduce the Landauer result. The orbital occupancies in the Landauer theory are defined by two chemical potentials (fixed by 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 figure 4(a). We see that our scheme leads to stable steady currents. We have tested convergence of our results with device size and the range of the density matrix from the device into the leads. For the device we allow only one atom to have the increased on-site energy with the remaining atoms being lead atoms. Results are summarized in figure 5. The current is generally below the Landauer value and is insensitive to the range of the density matrix from the device into the leads. As the reference density matrix is real, the imaginary part of the density matrix is suppressed by the damping term, thereby reducing the current in the leads. Damping dominates the electron dynamics when , where W is the bandwidth. In the present case .

Figure 5. The steady currents calculated as a function of the system size. Note that only one atom in the middle has the increased on-site energy of 1 eV. The bias is 0.1 V, the hopping integral is 1 eV and  eV.

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 Q_i, then it contributes a potential U_iQ_i (U_i  =  7 eV) to its own site and (where R_ij and U_i are in atomic units) to neighbouring sites. To match the potential at the ends of the device to that in the leads we include a linear external field in the device region of V (x)  =  a + bx. The results are shown in figure 4(b). It can be seen that, aside from a substantial 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, corresponding to a given temperature, and then perform molecular dynamics on it in one dimension. The position of the ion is allowed to evolve according to equation (4). Note that the environment atoms do not move, so no heat can be transported away.

To monitor ionic heating we follow the evolution of the kinetic energy with time (see figure 6). For a small bias (0.1 V) the ionic kinetic energy decays with time (cooling), while for a large bias (1.0 V) the kinetic energy increases (heating). The energy () 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 bias increases approximately tenfold the number of possible heating transitions, producing the observed changed behaviour.

This method is thus able to calculate some non-adiabatic effects of current flow. However, the amount of heating that we observe for a bias of 1 V is much less than we would expect from quantum perturbation theory [12]. The explanation for this, and a correction to remedy it, are the subject of ongoing work and will be presented in a paper currently in preparation.

The financial support of the IRC on nanotechnology for APH and of the Royal Society for DRB, and conversations with Tchavdar Todorov, Adrian Sutton and Marshall Stoneham are gratefully acknowledged.

References

[1]

Reed M A, Zhou C, Muller C J, Burgin T P and Tour J M 1997 Science 278 252 
CrossRef

[2]

Di Ventra M, Pantelides S T and Lang N D 2000 Appl. Phys. Lett. 76 3448 
CrossRef

[3]

Landauer R 1989 J. Phys.: Condens. Matter 1 8099 
IOPscience

[4]

Lang N 1995 Phys. Rev. B 52 5335 
CrossRef

[5]

Ness H and Fisher A J 1997 Phys. Rev. B 56 12469 
CrossRef

[6]

Sanvito S, Lambert C J, Jefferson J H and Bratkovsky A M 1999 Phys. Rev. B 59 11936 
CrossRef

[7]

Taylor J, Guo H and Wang J 2001 Phys. Rev. B 63 245407 
CrossRef

[8]

Brandbyge M, Mozos J L, Ordejon P, Taylor J and Stokbro K 2002 Phys. Rev. B 65 165401 
CrossRef

[9]

Todorov T N 2002 J. Phys.: Condens. Matter 14 3049 
IOPscience

[10]

Di Ventra M and Lang N D 2002 Phys. Rev. B 65 45402 
CrossRef

[11]

Montgomery M J, Todorov T N and Sutton A P 2002 J. Phys.: Condens. Matter 14 5377 
IOPscience

[12]

Todorov T N 1998 Phil. Mag. B 77 965 
CrossRef

[13]

Kohen D, Stillinger F H and Tully J C 1998 J. Chem. Phys. 109 4713 
CrossRef

[14]

Billing G D 1999 J. Chem. Phys. 110 5526 
CrossRef

[15]

Kapral R and Ciccotti G 1999 J. Chem. Phys. 110 8919 
CrossRef

[16]

Di Ventra M and Pantelides S T 2000 Phys. Rev. B 61 16207 
CrossRef

[17]

Goringe C M, Bowler D R and Hernández E 1997 Rep. Prog. Phys. 60 1447 
IOPscience

[18]

McWeeny R 1960 Rev. Mod. Phys. 32 335 
CrossRef

[19]

He L and Vanderbilt D 2001 Phys. Rev. Lett. 86 5341 
CrossRefPubMed

[20]

Goedecker S 1999 Rev. Mod. Phys. 71 1085 
CrossRef

[21]

Baer R and Neuhauser D 2003 Chem. Phys. Lett. 374 459 
CrossRef

[22]

Calvayrac F, Reinhard P G, Suraud E and Ullrich C A 2000 Phys. Rep. 337 493 
CrossRef

[23]

Horsfield A P and Bratkovsky A M 2002 J. Phys.: Condens. Matter 12 R1 
IOPscience

[24]

Todorov T N 2001 J. Phys.: Condens. Matter 13 10125 
IOPscience

[25]

Todorov T N 2003 private communication

[26]

Williams A R, Feibelman P J and Lang N D 1982 Phys. Rev. B 26 5433 
CrossRef

[27]

López Sancho M P, López Sancho J M and Rubio J 1984 J. Phys. F: Met. Phys. 14 1205 
IOPscience

[28]

Baer R and Neuhauser D 2003 Int. J. Quantum Chem. 91 524 
CrossRef