UFDC Home  myUFDC Home  Help 



Full Text  
TRANSPORT PROPERTIES AT NANO SCALES VIA FIRST PRINCIPLES STUDIES By CHUN ZHANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 Copyright 2004 by Chun Zhang To my parents. ACKNOWLEDGMENTS First, I would like to thank my advisor, Dr. HaiPing Cheng. She brought me into the field of computational physics and exposed me to its frontier. She continuously supported me throughout the past four years and most importantly, taught me that independent thinking is true heart of the research. Research is lonely and boring, but thanks to her, it can also be interesting. I would also like to thank Professors Arthur Hebard, Selman Hershfield, Adrian Roitberg and Samuel B. Trickey for serving on my supervisory committee. Dr. Xiaoguang Zhang, a senior scientist at Oak Ridge National Lab, is another important person to my Ph. D study. Not only his knowledge, but also his attitude toward the research will have a longterm impact to my career. My life at the University of Florida would have been black and white without my friends. They are Dr. Wuming Zhu, Dr. Maohua Du, Dr. Linlin Wang, Dr. Ping Jiang and Dr. Yao He. I would also like to extend special thanks to the following professors at Fudan University, Dr. Yuanxun Qiu, Dr. Jingguang Che and Dr. Jian Zi. Becoming their friend is one of the most important reasons that I want to be a physicist. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ................................................................................................. iv LIST OF FIGURES ......... ..................................... ....... vii A B ST R A C T .......... ..... ...................................................................................... x CHAPTER 1 O V E R V IEW ................................. .. .................................................... 1 2 FIRST PRINCIPLES APPROACH COMBINING NGFT AND DFT.....................4 2.1 M odel Hamiltonian and Landauertype Formula ..................................................4 2.2 The Combination of DFT and NGFT for An Open System ...................................8 2.2.1 SelfConsistent Procedure for Ground State of Closed System ...................8 2.2.2 Selfconsistent Procedure for Open System ........ ........... .................. 10 3 COHERENT ELECTRON TRANSPORT THROUGH AN AZOBENZENE MOLECULE: A LIGHTDRIVEN MOLECULAR SWITCH............... ...............12 3.1 Introduction ................................................. ........ .. ....................... 12 3.2 Tightbinding M odel and Selfenergies...................... ..... ................. 14 3.3 C onductance at Z ero B ias......................................................................... ... ... 18 3.4 The Sw itch under A Finite Bias....................................... ......................... 23 4 QUANTUM TUNNELING THROUGH MAGNETIC JUNCTIONS ......................30 4 .1 Introdu action ................ .............. ....... ...... ..................................... 30 4.2 Nonequilibrium Green's Function Formalism in Real Space ...........................34 4.3 Selfconsistency and Tunneling Current for A TMR Junction...........................37 4.4 Application to the Fe/FeO/MgO/Fe Tunneling Junction.............................. 43 4.5 Conclusion ..................... .................... ...........................61 5 PROBLEMS IN THE CURRENT NONEQUILIBRIUM APPROACH FOR OPEN SYSTEMS BASED ON DFT .............. .............. .... ............... 63 5.1 Introduction .................. ...... ............ ....................... ................ .....63 5.2 HohenbergKohn Theorem for Open System under Finite Bias ........................65 5.3 Exchange Energy for Uniform Electron Gas............. .. ........... ............... 75 5.4 Conclusion ..................................... ................................ .......... 82 6 SUMM ARY AND CONCLUSIONS..................................................................... 83 APPENDIX A EXPANSION OF LESSER GREEN"S FUNCTION ....................... ................85 B DETAILS OF NONEQUILIBRIUM SELFCONSISTENT PROCEDURE ...........88 L IST O F R E F E R E N C E S ......... .. ............... ................. ................................................90 BIOGRAPH ICAL SKETCH ...................................................... 94 LIST OF FIGURES Figure p 11. Partition of the open system connected to two metal electrodes ............................2 31. Physical models of an azobenzene switch in contact with two gold leads: trans isomer (upper) and cis isomer (lower). ........................................ ............... 14 32. Tightbinding model of an open system connected to two semiinfinite leads. Indices i, j denote nodes of device region and a, b are the nearest node of left and right lead. Only the nearestneighbor interaction is considered. .....................................15 33. A perfect lead that consists of two semiinfinite leads and a hopping term. Indices m and i belong to the left semiinfinite lead and i+1, j belong to the right semiinfinite lead ................................................................................. 16 34. Density of states (DOS) projected onto the switch region (upper panel) and transmission as a function of energy (lower panel). The solid lines represent the trans configuration and dotted lines the cis. Intense peaks in the DOS are truncated to im prove readability. ......................................... ......................... 19 35. Local DOS of the trans and cis configurations decomposed into contributions on the Au layers (lower panel), S atoms (middle panel), and the azobenzene molecule (upper panel). The solid lines represent the trans configuration and dotted lines cis.20 36. Local DOS at the Fermi energy projected on seven sites in the switch. From 17, or from left to right, Au lead, S atom, benzene ring, doublebonded N2, another benzene ring, S atom and Au lead ............................... .....................22 37. Nonequilibrium selfconsistent procedure. The nonequilibrium Green's function is calculated according to equation (211). ....................................... ............... 24 38. The difference of electronstatic potential between bias one volt and bias zero for cis configuration. The xz plane is parallel to the N2 bond, cutting through the m idpoint of A u unit cell. ................................................................. .....................25 39. The transmission probability as a function of singleparticle energy and the external bias voltage for cis configuration. ........................................ ......................... 26 310. Current as a function of bias voltage for both trans and cis configurations. The red line represents the current of trans and the blue line represents the current of cis...27 311. The differential conductance as a function of bias voltage for both trans and cis configurations............................................................................................. 28 312. The transmission as a function of energy at 0.64, 0.72 and 0.8 Volts for cis configuration ...................................................... ................. 2 8 41. Shifts in the density of states in each atomic layer at a bias voltage of 0.54 V..........44 42. Derivative of the capacitive charge with respect to the bias voltage. The charge on the right electrode includes that on the oxygen atom in the interface MgO layer. 'Net charge' includes that on both the electrode and the dielectric charges on the MgO layers on the same side of the midpoint of the barrier...............................46 43. Derivative of the moment with respect to the bias voltage. The labels 'left' and 'right' indicate that the layers to the left (right) of the midpoint of the barrier are summed. The contribution to the sum from the nonbulk layers of the iron electrodes is taken to be the difference between the moment of that layer and the bulk moment in the electrode on the sam e side ....................................... ................. ......... 47 44. Zero bias transmission probability as a function of k/, (a) Majority spin channel for parallel alignment of the moments; (b) The dominant spin channel for antiparallel alignment of the moment......................................... ...... ............... 50 45. Current density as a function of bias voltage for majority and minority spin channels for the case of parallel alignment of the moments between electrodes ..................51 46. Current density as a function of bias voltage for majority and minority spin channels for the case of antiparallel alignment of the moments between electrodes.............51 47. Transmission probability at energy E = /, as a function of k// along the line of k = ky for the majority spin channel for parallel alignment of the moments at different bias voltages. ..................... ................ .............................53 48. Transmission probability at energy E = /, as a function of k// along the line of k = ky for the minority spin channel for antiparallel alignment of the moments at different bias voltages. ..................... ................ .............................55 49. Integrated transmission probability as a function of k// along the line of kx = k, for the majority spin channel for parallel alignment of the moments at different bias v o lta g e s ....................................................................... 5 6 410. Integrated transmission probability as a function of k// along the line of k = ky for the majority spin channel for parallel alignment of the moments at different bias v o lta g e s ....................................................................... 5 8 411. Integrated transmission probability as a function of k// for the majority spin channel for parallel alignment of the moments at bias 0.378 V. ................. ...............59 412. Integrated transmission probability as a function of k// for the minority spin channel for antiparallel alignment of the moments at bias 0.378 V. ................................59 413. Magnetoresistance ratio as a function of bias voltage. Point A denotes the MR ratio at zero bias by removing all contribution from interface resonance states .............61 51. Exchange energy at different bias voltages for different p (Electrons/m3) .............78 52. Split of the exchange energy: (a) the exchange between nonequilibrium electrons and the exchange between nonequilibrium and equilibrium electrons and (b) the exchange between equilibrium electrons. ..................................... ............... 80 53. Exchange energy as the functional of p/ p for different electron densities ..........81 A1. Expansion of lesser Green's function. The time of the perturbation can be either on + branch or on branch...................... .. .. ......... .. ........................ ............... 86 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy TRANSPORT PROPERTIES AT NANO SCALES VIA FIRST PRINCIPLES STUDIES By Chun Zhang December 2004 Chair: HaiPing Cheng Major Department: Physics There are two main difficulties for the first principles study of transport properties at the nano scale. The first is that manybody interactions need to be taken into account for the infinite system without periodic boundary conditions. The other is that the system is usually in a nonequilibrium state. Both of these two difficulties are beyond the ability of conventional first principles methods to reconcile. Recently, a new first principles approach which combines the Nonequilibrium Green's Functions Technique (NGFT) and the Density Functional Theory (DFT) was proposed. DFT has been proved to be successful in molecular and solid state physics. Currently used DFT approximations can take into account 'most' manybody effects and NGFT naturally includes the non equilibrium effects. The new approach uses NGFT to treat the nonperiodic boundary conditions and DFT to treat many body interactions. This approach has been successfully used in molecular electronics. The thesis is organized in the following way. First: we introduce the main ideas of combining NGFT and DFT, and then apply this method to a lightdriven molecular switch. The switch, made of a single molecule, is one of the most important elements of nanoelectronics. However, most proposed molecular switches are driven either by an external bias voltage or by STM manipulation, neither of which is ideal for nanoscale circuits. The switch we designed has a high onoff conductance ratio and more importantly, can be driven by photons. In following chapters, we generalize the method to the spindependent case and apply it to a magnetic layered structure. We implemented the method within the framework of the Layer KorringaKohnRostoker (LKKR) approach, which is particularly welladapted to the layered structure and found a biasenhanced tunneling magnetoresistance (TMR) for the Fe/FeO/MgO/Fe junction. Our results are important not only for application, but also for understanding of the voltagedependence of TMR for layered structures. The experimental studies show that the bias voltage usually kills the TMR of amorphous magnetic tunneling junctions. Our study shows that for an impurityfree layered structure, a different behavior of TMR may occur. CHAPTER 1 OVERVIEW As the trend of miniaturization continues, we eventually enter the regime of nano electronics, the technology of building electronic devices on the basis of a single molecule. Experiments have shown the great potential of modem technologies, such as STM and AFM, to make molecular transistors and switches etc. [15]. For theoretical studies, the circuit in which the molecule is embedded provides great difficulties. First, the molecule is connected to the circuit through two open boundaries. It is crucial to describe the open boundaries in a way that the interaction between the molecule and the complicated circuit is taken into account. Second, the circuit is almost always far away from the equilibrium state and in most cases, varying with time. The timedependent process is not the focus of this thesis. We assume that the whole circuit is in a steady state. Our goal is to find the steadystate solution for an open system far away from equilibrium. An open system connected to two metal electrodes usually has geometry such as shown in Fig. 1. The whole system can be divided into three parts, the left lead, the right lead and the device region. The device region includes the molecule and enough portions of the two metal electrodes so that the electrons in the two leads are not affected by interactions in the device region because of screening. A finite bias is applied to the system. The goal of the research is to relate the steadystate current flowing through the device region to the external bias voltage. Traditionally, the electron transport phenomena through such an open system are studied by the Boltzman equation in conjunction with the effectivemass assumption. The Boltzman equation is based on a probabilistic description for the evolution of a singleparticle distribution f(ri, p, t) which is a function of classical variables r and p However, when the device region is at nanoscale, the phase coherence of charge carriers is usually maintained in the whole device region so that quantum effects dominate. For this kind of coherent quantum transport phenomena, the first principles approach which combines the nonequilibrium Green's functions technique (NGFT) and Density Functional Theory (DFT) has proven to be successful [622]. The method can be summarized in two steps. First, a non equilibrium selfconsistent procedure based on DFT and NGFT is used to determine the singleparticle effective potential. Second, a Landauertype formula is used to calculate the current, which relates the current to the integral of the singleparticle transmission probability over energies. Current I Left lead Device Region Right lead Figure 11. Partition of the open system connected to two metal electrodes. The approach combining NGFT and DFT is based on the singleparticle picture. The tunneling process is explained within this framework as follows. The Bloch electrons in two leads tunnel through the molecular orbital in the device region. The success of this picture can be best seen in the work recently published by a Denmark group [21], where the current flowing through a broken Au nanowire is calculated. At small bias voltages, the theory agrees well with the experiment. The approach can also be equivalently built on a scattering wave formulation [8 11]. The scattering wave approach has so far been implemented only with jellium electrodes, due to the difficulties associated with computing and storing Bloch wave functions of realistic leads. In Chapter 2, we introduce the details of the twostep approach. In Chapter 3, we apply the method to the Au/Azobenzene/Au junction. Our calculation shows that the azobenzene molecule is a good candidate for the lightdriven molecular switch. In Chapter 4, we generalize the method to spindependent tunneling junctions with application to a layered magnetic Fe/FeO/MgO/Fe junction. The calculation suggests a biasenhanced tunneling magnetoresistance (TMR), which is important for both application and the understanding of TMR through a layered structure. In Chapter 5, we discuss the problems that exist in the current twostep approach. For the first step, in which the effective singleparticle potential is determined, DFT is the biggest problem. The HohenbergKohn theorem, the basis of DFT, is correct only for the ground state [23]. The use of DFT for the open system far from equilibrium state is problematic. We discuss this problem in detail and present a possible way to introduce nonequilibrium effects into DFT. For the second step, in which a Landauertype formula is used to calculate the current, the problem is that the KohnSham orbitals are used to calculate the current. Such orbitals are not genuine physical quantities, which causes trouble in choosing basis sets and functionals. CHAPTER 2 FIRST PRINCIPLES APPROACH COMBINING NGFT AND DFT 2.1 Model Hamiltonian and Landauertype Formula As we mentioned in the previous chapter, the system of interest can be divided into three parts, the left lead, the device region and the right lead. To relate the current flowing through the device region to the external bias is the goal of our research. The whole system is infinitely large, but without periodic boundary conditions, which provides a big challenge to the theoretical study. However, if one assumes that the two leads are noninteracting, an exact formula which relates the current to local properties of the device region can be derived [24]. The derivation is outlined in following paragraphs. The starting point is the Hamiltonian H = H + R + Hce, + HoC (21) where HL HR and HCe,, denote the Hamiltonian for the left lead, the right lead, and the device region respectively. The last term HCo, denotes the coupling between leads and the device region. Since the two leads are noninteracting, we have HLR = kak,a k,aeL,R en = nt({d, {d) (22) Hfou, = IZ[Vk,cT ldn + h.c. k,aEL,R n where indices k, a denote the channel of left and right leads (Bloch wave vectors and spin) and n denotes the singleparticle basis of the device region. The key issue is how to treat the interaction of the device region. The expression for the current can be derived by taking the time derivative of the total number of electrons in the left lead J = e(N )ie([H, NL (23) where NL = 'cck, and JL is the current from the left lead into the device region. k,cxL Inserting the Hamiltonian into the expression for the current, we have JL =ie Vkn (ca d,) Vk~ (dc k,aGL (24) e knka Vka,nG k,aGL n where the current is expressed in terms of two nonequilibrium 'lesser' Green's functions. Now, the problem becomes how to calculate the 'lesser' Green's functions. This can be achieved by a perturbation procedure. At t = o, the two leads and the device region are well separated, which means that there is no the coupling term in the Hamiltonian. The three parts maintains their own Fermi energies. Without loss of generality, we assume that the Fermi energy of the right lead /Ru is higher than the Fermi energy of left lead/L,. As we gradually turn on the coupling from t = c tot = 0, the system will achieve a steady state at t = 0 with a current J from left to right. Taking the coupling term in the Hamiltonian to be the perturbation, the 'lesser' Green's functions can be expanded as Gn,ka ka,m[G nmka,ka nmk+G ,ka n (25) ka,,n = Vka,m ka,kaG n k+ gakaM n m where m runs over all basis states in the device region g denotes the Green's function for noninteracting leads and the upper indices r, a denote retarded and advanced Green's functions respectively. The details of the expansion can be found in Appendix A. Introducing the equation (25) into equation (24) and doing a Fourier transform gives J, = 2e d Re{ VknVk m k ()Gm (0) + gK~ ,ka ()Gn( (o)]} (26) f2;r kaeL n,m where the summation runs over all channels of the left lead. Since the lead is non interacting, the Green's functions g g" can be calculated as gka,ka (0) = 27if(Eka )5(O) ka) 1 (27) gka,ka ka i where f(ska) is the FermiDirac distribution function. Insertion of the equation (27) into equation (26) gives the current as J = ie Tr{L [fL (Gr G) + G]} (28) f 2;r where the linewidth function FL is defined as [F (k )]mn = 2r ( pa (c )V. (,k )Vl (k) (29) acL in which the term pa is the density of states in the lead. A similar formula can be derived for J, the current from the right lead into the device region, so the total current is given 1 byJ (JL JR),whichis 2 = ie dTr[FL FR ]G +[fLFL fRFR][Gr G]} (210) 2 f 2;r where the current is related to the local properties of the device region. In equation (210), the Green's functions in the device region are solvable for only a few models. The simplest one among them is the free lectron model, which means that the device region is also noninteracting. For the freeelectron model, we have [15, 20, 25] G' = Gr''G Y< = i[F/f, +FRR (211) and then the current can be written as J= elde[f fR ]Tr{GTRGTr } (212) which is a Landauertype formula. Since the interaction in both leads and the device region usually plays an important role in the tunneling process at nano scales, equation (212), which is based on the free electron model, seems to be uninteresting, unless it is combined with a mean field theory. DFT is one kind of mean field theory. It is based on the singleparticle picture. Within presentday approximations, it is able to take into account 'most' manybody effects. Moreover, due to the success of DFT in the electronic structure calculation of molecules, it is natural to combine DFT and the Landauertype formula to study molecular electronics. 2.2 The Combination of DFT and NGFT for An Open System 2.2.1 SelfConsistent Procedure for Ground State of Closed System DFT has been known to be powerful in the electronic structure calculation for both molecular and solid state systems [2630]. The main ideas of DFT are outlined in the following paragraphs. The basis of DFT is the HohenbergKohn theorem [23]. This theorem shows that the ground state of a closed system is uniquely determined by the electron density of the system, which means that average value of any physical quantity for the ground state can be expressed as a functional of the electron density. As an example, the total energy can be written as E[p] = F[p] + p()Vext (i)d (213) where the term Vex, represents the entire external potential (including the ionelectron interaction) and the first term F[p] is a universal functional related to the electron electron interaction and the kinetic energy of electrons. The theorem also shows that the true groundstate electron density minimizes the energy functional E[p]. The only thing used in the proof of HohenbergKohn Theorem is the minimum energy principle of ground state. On the basis of the HohenbergKohn theorem, a selfconsistent procedure based on the singleparticle picture can be constructed. Assume that there is a noninteracting system which yields the same electron density as the interacting system of interest, then the electron density can be written in terms of the singleparticle functions p(/) = ~ (F) 2 (214) 1 which are usually called KohnSham orbitals. The kinetic energy of the noninteracting system can be calculated in terms of the KohnSham orbitals (Hartree atomic units) V2 T/[pT]=  ). (215) Then the exchangecorrelation energy can be defined as Exc [p] = F[p] T,[p] EH [p] (216) where the term EH [p] is the classical Coulomb energy. The exact form of Ex[p] is unknown. However many kinds of exchangecorrelation functionals, most based on the local density approximation, have been proposed and successfully used. By minimizing the energy functional, we get the KohnSham equations V2 + VH [p](r) + Ve, () + V [P](ir), (i)= g () (215) where the first term VH is the wellknown classical Coulomb potential and the term Vxe is the exchangecorrelation potential defined by [ The effective KohnSham 9p potential is defined as VeY [p]() = VH [p](r) + Ve, (r) + VX [p]() (216) Now, a selfconsistent procedure (KohnSham procedure) can be completed as follows. Given the electron density p, the effective KohnSham potential Vff [p] can be calculated according to equation (216). Knowing VJ. [p], we can get the KohnSham orbitals V, by solving equation (215) and then the new electron density can be calculated according to equation (214). Repeat the procedure until the selfconsistency is achieved. The HohenbergKohn theorem guarantees that the KohnSham procedure yields the true ground state electron density (the V representability problem for the electron density is not discussed in this thesis). Thus an interacting system can be solved by a selfconsistent procedure based on the singleparticle picture. 2.2.2 Selfconsistent Procedure for Open System For an open system with the geometry shown in Fig. 11, there are two kinds of eigenstates, scattering states which are extended in all space and bound states which are localized in the device region. If the bias voltage is zero, DFT is still applicable due to the fact that the system is indeed in a ground state (the temperature is set to zero). In practice, since we cannot treat infinitely large system, certain approximations have to be made. The major approximation is that, due to screening effects, the electron density and effective KohnSham potential in the two leads are not affected by the device region and thus can be set to the bulk values. The bulk is also infinitely large, but with the help of the Bloch theorem, we only need to do the DFT calculation in one unit cell. After fixing the effective KohnSham potential of the two leads, the remaining problem is to treat the device region, which is finite, with two fixed boundaries. If the bias voltage is not zero, the use of DFT becomes problematic. In Chapter 5, we will discuss this problem in detail, but here, we just assume that the same approach can be generalized to the finite bias case. The goal is to calculate the electron density selfconsistently in the device region with two fixed, open boundaries. The Green's functions technique described below is an efficient way to do this. The electron density in the device region is directly related to the 'lesser' Green's function by p = G< (c)d (217) 2n i For a singleparticle Hamiltonian, the 'lesser' Green's function can be calculated by the equation (211), where the retarded and advanced Green's functions are related to the effective potential of the device region and two open boundaries as Gr(a) = 1 (218) E C y r(a) r In equation (218), the term HC,, denotes the effective singleparticle Hamiltonian and the two terms EL ,R are selfenergies for the left and right lead respectively. Two self energy terms represents effects of open boundaries. We will see an example of calculating these two terms in Chapter 3. As a conclusion, the selfconsistent procedure for an open system can be summarized as follows. First, the two leads can be treated by a separate bulk calculation which provides two open boundaries for the device region. This only needs to be done once. Second, given the electron density p in the device region, the effective Kohn Sham potential can be calculated. Knowing the effective potential, the 'lesser' Green's function can be computed and then a new electron density is given by equation (217). After selfconsistency is achieved, the Landauertype formula (212) can be used to calculate the current in the system. The nonequilibrium selfconsistent procedure mentioned above is implemented in a Green's functions representation which employs a multiband tightbinding model in Chapter 3 and is implemented in a wave function representation which includes spin in Chapter 4. CHAPTER 3 COHERENT ELECTRON TRANSPORT THROUGH AN AZOBENZENE MOLECULE: A LIGHTDRIVEN MOLECULAR SWITCH 3.1 Introduction Electronic devices soon will be constructed of molecules or have features of molecular size. The results of recent experimental [15] and theoretical studies [620] suggest a brilliant future for molecular electronics. Among various efforts and activities, several schemes have been proposed to design and construct molecular switches that have two or more distinct states with vastly different conductance. Switching between on and off states can be performed by applying an external bias or by using a scanning tunneling microscope tip to manipulate the system [3135]. These methods are not ideal. The former approach may interfere greatly with the function of a nanosize circuit, while the latter one imposes severe limitations in device applications. In addition, these methods are relatively slow. Alternative methods, such as those based on fast, lightdriven processes, are therefore highly desirable. The photosensitive azobenzene molecule discussed here is one such candidate for an optoelectronic device. The azobenzene molecule has two stable conformations [3638], trans and cis. The energy of the trans structure is 0.6 eV lower than the energy of the cis structure [39], and they are separated by a barrier of approximately 1.6 eV (measured from trans to cis) [40]. The trans structure is 1.98 A longer than the cis. The molecule can switch reversibly from one structure to another under photoexcitation, with the structural change occurring on an electronically excited state [38, 41]. Previous work has shown that a nanoscale optomechanical device based on the conformation change of the azobenzene polymer can be designed. Experimentally [41], it is known that light with a wavelength of 365 nm transforms azobenzene units from the trans to the cis structure. As a result, the polymer contracts, thus converting the energy of electromagnetic radiation to mechanical work. A second pulse of light with a wavelength of 420 nm reverses the cis to the trans configuration. In this chapter, we report investigations of transport in the molecular wire, AuS azobenzeneSAu. It consists of a single azobenzene molecule, which has been functionalized by replacing a hydrogen atom with a CH2S group to provide a good contact with gold, and two semiinfinite gold leads. Note that the name of the modified molecule is technically bis(4methanethiolpheyl) diazene or BMTPD. Since one can add more azobenzene units to elongate the polymer, we simply follow Refs. [38, 39] and refer to the molecular system as 'azobenzene'. As shown in Fig. 31 and consistent with previous chapters, the entire system is divided into three regions: left lead, switch, and right lead. The switch includes two layers of gold from each lead. The structures of the leads are chosen to be the same as in Ref. [42], which corresponds to a gold (111) crystal lattice. The electronic structure of the switch region is determined by treating it as a cluster. Full quantum calculations, based on density functional theory (DFT), are used to optimize the structure of a cluster consisting of the molecule, the sulfur atoms, and two rigid layers of gold. The optimization is performed with two different choices of the interlead distance, d. For the trans configuration, the calculated optimal d is 18.16 A and, for the cis configuration, 16.23 A. In subsequent sections, we study the transport properties for both equilibrium and nonequilibrium case. In both cases, the localized basis set and the tightbinding model is employed. k=365 nm X=420 nm Figure 31. Physical models of an azobenzene switch in contact with two gold leads: trans isomer (upper) and cis isomer (lower). 3.2 Tightbinding Model and Selfenergies As mentioned in Chapter 2, the key issue is to evaluate the self energies of two semiinifinite leads, 2L and XR Once these two terms are computed, we can construct the retarded and advanced Green's functions G' and G" of the device region in the presence of two leads. Then the lesser Green's function can be calculated according to equation (2 11). The selfenergy term represents the effects of the semiinfinite lead on the device region. Usually, a truncation of the Hamiltonian in real space is needed since we can not treat an infinitely large system. So, the tightbinding model is a natural choice. Fig. 32 sketches the model, where indices i, j denote nodes of the device region and a, b are the nearest node of the left and right lead respectively. Note that hopping terms h, h' are operators [42]. h: a i J b Figure 32. Tightbinding model of an open system connected to two semiinfinite leads. Indices i, j denote nodes of device region and a, b are the nearest node of left and right lead. Only the nearestneighbor interaction is considered. Within the framework of tight binding, the Green's function G, can be written as G, = G ,c" + GC"h g hG, + GCh'gbh 'G, (31) where the two hopping terms are taken as perturbations and the Dyson equation is used. The term Gc"e represents the Green's function of the device region without two leads and the term g denotes the unperturbed Green's function of the leads. Then, self energies can be defined as S= h gh (32) R = h gbbh where, EL, 2R are selfenergies of the left and right lead respectively. Two hopping terms are projections of the effective KohnSham Hamiltonian to a localized basis of each node (hopping terms are matrices also). The only problem remaining is two socalled surface Green's functions ga and gbb Since a perfect lead with periodic boundary conditions consists of two semiinfinite leads and a hoping term within the framework of tightbinding as shown in Fig. 33, it is possible to construct the Green's functions of semiinfinite leads on the basis of Green's functions of a perfect lead by a perturbation procedure. The whole procedure is described in the following. h m i i+1 j Figure 33. A perfect lead that consists of two semiinfinite leads and a hopping term. Indices m and i belong to the left semiinfinite lead and i+1, j belong to the right semiinfinite lead. We take the hopping term h to be the perturbation, then the node i is the surface node of the left lead and the node i+1 is the surface node of the right lead. G is used to denote the Green's function of the perfect lead and Go, the Green's function of two semi infinite leads. The Dyson's equation reads G = G, +(i GVGj) = GOVGj (33) where p, q run over the whole system. Since V only contains two terms Vii+ and Vi+,i . P equals to i or i+1. Because Gii+l=0 (without perturbation, an electron can not propagate from i to i+1), we have G, = GyV,,+G,G+, = glhG,+, (34) Similarly, G,,, = G y G, = g+ ,zh G, ,. (35) In order to solve g from equation (34) and (35), we have to know relations between Gi and Gi+,j. This can be found by solving the Green's functions of the perfect lead. For a perfect lead, the tightbinding Hamiltonian is H = hoa,, + C haa,, + h a+, am, (36) m m where h,,h are N x N matrices, N being the number of orbitals within one unit cell of the lead. The eigenvalue problem for Bloch states takes the form [ h, (h E)ekd (h E*)ekd (k, e) = 0, (37) where k is the Bloch vector and d is the lattice constant. The energy term e contains an infinitesimal imaginary part, E' = E+ ir with r>> 0. Equation (37) can be written as Az h h z 2=, =0, (38) where, M+ = E+ ho and h, = h E It is to be noted that not all eigenvalues A, of equation (38) for a given energy represent the Bloch state. Because equation (38) is quadratic in A and has a dimension of N, the number of eigenvalues for a given energy is 2N. It is shown in reference [22] that, when the energy is not real, half of those eigenvalues are outside the unitary circle in the complex plane (we denote them as A> ) and another half are inside the unitary circle (we denote them as A< ). Using the superscripts + to denote eigenvalues and eigenvectors for e+ and E respectively, we can introduce two Nx N matrices, a and f, a = UA(U) U (39) 8+ = U+A* (U) ), where U+ is a matrix, whose i th column is the eigenvector u,, and A is a diagonal matrix, whose i th diagonal element is the eigenvalue A2 The matrices, a and f, satisfy the matrix equation Ma h (h)+(a' )2 = (310) MA/1 + (h) h (3+ )2 0 The retarded Green's function is defined by the matrix equation (E H)G+ = I. We have (i( H)G j)=8 M+G, h G,, (h )+ G+, = 0. (for i < j). (311) For a perfect lead, the solution of this equation is done via recursion relations of Green's functions, G = AG+, andGi, = AG G,. This gives the equation MA h (h ) A2 = 0. (312) Comparing with equation (310), we see that, for i < j, A = a Introducing G+ = AG+, into equation (35), we have g, = a (h) (313) Similarly, fori > j, we have g ,+, = p+ (h+ ) (314) 3.3 Conductance at Zero Bias In this section, we calculate the conductance of the two optimized structures. In this study, we focus on the linearresponse regime. Within this limit, the entire system is in equilibrium and there is a common Fermi energy across the whole system. Because of charge screening by electrons in the metal, the two leads, as well as the switch, are charge neutral. The conductance is calculated by the Landauer formula [44], 2e2 ar= T(E ) (39) h where T(EF) is the transmission function at the Fermi energy. The effective potentials are calculated by the computational chemistry code NWCHEM [28]. In the calculation, we use Crenbs effective core potential (ECP) Gaussian basis [28] for Au atoms, STO3G Gaussian basis [28] for all other atoms (C, S, N, H). The pbe0 exchangecorrelation functionals [26] are used. For the wire system depicted in Fig. 31, the system is, in principle, infinite, and periodic boundary conditions are not appropriate. However, due to screening, the KohnSham potentials in the two semiinfinite leads are not affected by the switch. As a result, the potentials can be computed in separate calculations of a one dimensional gold wire with periodic boundary conditions [20, 22, 42]. 400 U a Trans b b > 200 a SE, 08 S Cis E \ Trans S04 00, 024 022 020 018 Energy (Hartree) Figure 34. Density of states (DOS) projected onto the switch region (upper panel) and transmission as a function of energy (lower panel). The solid lines represent the trans configuration and dotted lines the cis. Intense peaks in the DOS are truncated to improve readability. With Green's function techniques, the effect of the two leads on the switch can be absorbed into selfenergy terms, which are calculated using the same method as mentioned in previous section. The transmission function T(EF) is evaluated by the Caroli formula [25] T= Tr[FLGrFRGa] (310) which is the kernel of the integral in equation (212). The density of states (DOS) in the switch region can be calculated from the trace of the retarded Green's function, which includes both scattering states and localized states. The DOS can then be projected onto any desired subsystem. Note that the subsystem can be chosen as a single atom or a group of atoms. The peaks in the DOS as a function of energy are often used as indications of the molecular orbital and, with the combination of the transmission as a function of energy, useful information can be extracted as to how any given molecular orbital contributes to the conductance. Azobenzene molecule Trans b b' 200 Ef Cis S 0 o Trans > 200 Sulpher atoms a SCis a a' Trans b' 200 Au layers Cis 024 022 020 018 Energy (Hartree) Figure 35. Local DOS of the trans and cis configurations decomposed into contributions on the Au layers (lower panel), S atoms (middle panel), and the azobenzene molecule (upper panel). The solid lines represent the trans configuration and dotted lines cis. The calculated transmission as a function of energy and the DOS for the trans and the cis structures, with optimal interlead distance d, are shown in Fig. 34. Unlike the DOS of an isolated molecule, Fig. 34 displays a series of peaks superimposed on a broad background. These features originate from a combination of localized states from the molecule and a nearly continuous energy band of the two leads. At some energies, the DOS shows peaks, but the total transmission is very small. Notice, for example, peaks a and b, for the trans isomer, and peaks a' and b' for the cis isomer. These peaks are indications of localized states, which can be seen more clearly in the local density of states. As depicted in Fig. 35, the DOS of the switch region is decomposed into three subsystems: the two Au layers, the two S atoms, and the azobenzene molecule. For the trans isomer, peak a is mainly localized on the contacts, namely the Au layers and S atoms, while peak b is localized mainly on the molecule. For the cis isomer, peak a' is localized mainly in the Au layers, while peak b' has high density on the molecule. The DOS on the S atoms is very small in the cis isomer as compared to the same peak in the trans isomer. These comparisons indicate that the trans and cis configurations have quite different contacts with the leads. Furthermore, Fig. 35 indicates clearly that the broad background in the DOS derives primarily from the Au layers. Figure 35 shows that the Fermi energy of the system lies between the highest occupied molecule orbital (HOMO) and the lowest unoccupied molecule orbital (LUMO) of the azobenzene molecule, slightly closer in energy to the HOMO. Compared to the isolated molecule, the peaks in the local DOS projected on the molecule indicate that the molecular orbitals are perturbed by the presence of the leads and contribute to scattering states near the Fermi energy. Since we focus here on the zerobias case, the conductance is determined completely by transmission at the Fermi energy, which is dominated by the tail of the electron distribution tunneling through the perturbed HOMO. Note that the local DOS of both the Au layers and the S atoms peak near the trans HOMO energy, but that the magnitudes of these peaks are minimal for the transmission function in Fig. 34. For the trans isomer, tunneling through the perturbed HOMO leads to a very broad and intense peak in the transmission function, but a sharp and smaller peak for the cis isomer (Fig.34, lower panel). As a result, tunneling at the Fermi energy is more rapid through the trans structure than through the cis structure, although the DOS of the cis isomer is higher than that of the trans isomer. The calculated zerobias conductance for the trans structure is 0.21 x 104 41 which indicates a moderately good conductor. In contrast, the conductance for the cis structure is nearly 2 orders of magnitude lower. 30 1: Left Au layer 2: Left S atom 25 3: Left benzene ring S4: N=N 20 Cis 5: Right benzene ring S6: Right S atom >, 7: Right Au layer S 15 10 ) \ 1 5   0 Trans U \ 0 0 1 2 3 4 5 6 7 Atom groups Figure 36. Local DOS at the Fermi energy projected on seven sites in the switch. From 17, or from left to right, Au lead, S atom, benzene ring, doublebonded N2, another benzene ring, S atom and Au lead. To understand the tunneling process at the Fermi energy, Fig. 36 depicts the local DOS at the Fermi energy, Fig. 36 depicts the local DOS at the Fermi energy for both the trans and the cis isomers. In this figure, we decompose the DOS into seven components: the two left Au layers, the left S atom, the left benzene ring with a CH2 group, the doublebonded N2, the right benzene ring with a CH2 group, the right S atom, and the two right Au layers. As can be seen in the figure, the local DOS for the two molecular conformations differs significantly. For the cis isomer, the local DOS is much higher on the Au layers and the two benzene rings than on the S atoms and the N double bond. The nearly zero DOS on the right S atom and the oscillation of the DOS along different sites suggest that the scattering states are localized primarily on the benzene rings and the gold leads. Reflection of the scattering wave occurs at sites of low DOS (on S and N). Both localization and reflection contribute significantly to the resistance of the device. Compared to the cis isomer, the local DOS for the trans isomer displays a smoothly varying curve, suggesting that localization and reflection at the N double bond and the two S atoms contribute much less to resistance than in the cis configuration. In the azobenzene/gold system, the extended scattering states allow facile electron tunneling through the molecule. This result can also be interpreted from the viewpoint of molecular structure. In the cis configuration, due to the different orientations of the two benzene rings, orbitals on different rings have different configurations. Consequently, when electrons tunnel through the azo unit of the cis configuration, they are scattered more than they are in the trans configuration. A more illuminating way to understand the tunneling process is to decompose the conductance into conduction channels of two leads, which will be done in future work. The conductance of the cis isomer is significantly smaller than that of the trans. Since the azobenzene molecule can easily change between trans and cis configurations via excitation with light, we suggest that this system is a viable candidate as a light driven molecular switch. The switch will work at room temperature due to the high thermal stability of both the trans and the cis configurations of the azobenzene molecule. 3.4 The Switch under a Finite Bias When a finite bias is applied, the whole system is out of equilibrium. In this section, we apply the nonequilibrium selfconsistent procedure discussed in Chapter 2 to the Au/S/azobenzene/S/Au junction and relate the current through the switch to the external bias voltage. In the study, we assume that the geometry of the device region is the same as that of zerobias case. We neglect the effects of the socalled currentinduced force by making this assumption. 1 G 2; i Vff = V[p] Figure 37. Nonequilibrium selfconsistent procedure. The nonequilibrium Green's function is calculated according to equation (211). The selfconsistent procedure is depicted in Fig. 37. Since we are working within the framework of DFT, the most important physical quantity is the electron density. Given the electron density p, the effective potential Ve, can be computed as a functional of the density. Once we know the effective potential, the nonequilibrium Green's function G' can be calculated according to equation (211), where the retarded and advanced Green's functions are calculated in the same way as in the previous section. The details of the selfconsistent procedure can be found in Appendix B. One of the major accomplishments of the nonequilibrium first principles calculation is that the electrostatic potential due to the external bias can be computed self consistently. Fig. 38 shows the difference of the electrostatic potential between bias one volt and zero bias for the cis configuration. In the figure, the y axis is along the current direction and the xz plane is parallel to N double bond of cis configuration, cutting through the midpoint of the Au unit cell. The drop of the electrostatic potential across the switch area can be clearly seen. Cis under 1 Volts 0 0.5 . '0.0 01 S(Angstros 10 Figure 38. The difference of electronstatic potential between bias one volt and bias zero for cis configuration. The xz plane is parallel to the N2 bond, cutting through the midpoint of Au unit cell. The current through the system is calculated by equation (212), where the current is expressed as the integral of the singleparticle transmission probability over energy. The effects of external bias voltage on the current through the switch enter in two ways. First, as the bias voltage increases, the upper and lower limit of the integral changes so that more channels enter the integral. Second, the bias voltage changes the electrostatic potential, thus changes the transmission probability in turn. In Figure 39, we show the transmission probability as a function of both energy and bias voltage, where the transmission is expressed as T(E, V) = Tr(FL (E, V)G' (E, V)F (E, V)G' (E, V)). (313) Transmission as a function of energy and voltage for cis 1.0 HOMO c / LUMO 0.0Ie 7 0\ \ nergyV) 4 00 ,9 Figure 39. The transmission probability as a function of singleparticle energy and the external bias voltage for cis configuration. There are two effects of the external bias on the transmission probability through a particular molecular orbital. First, as the external bias changes, the effective KohnSham potential in the device region changes due to the change of electrostatic potential, thus the 'molecular orbital' changes. On the other hand, the external bias shifts the Fermi energies of the two leads in opposite directions so that channels in both electrodes corresponding to the same molecular orbital change. Both of these two effects are taken into account in the first principles study. IV curves through the trans and cis configuration are shown in Fig. 310. The differential conductance, defined as the derivative of the current with respect to the bias voltage, is shown in Fig. 311. One of the important results is that at zero bias, the conductance of trans configuration is two orders higher than the cis, which agrees with the equilibrium calculation. The current through the trans configuration is much higher than the cis for all bias voltages, which shows the high stability of the switch under a finite bias. 25 20 a15 S0 Trans SCis 5 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Voltage (Volts) Figure 310. Current as a function of bias voltage for both trans and cis configurations. The red line represents the current of trans and the blue line represents the current of cis. An interesting feature is shown in Fig. 310. For the cis configuration, when the bias increases from 0.7 Volts to 0.8 Volts, the current drops, which gives a negative differential conductance as shown in Fig. 311. This can be explained by consideration of Fig. 312, where the transmission as a function of energy at 0.64, 0.72 and 0.8 Volts are shown. As we mentioned before, as the bias voltage increases, more channels enter the integral of equation (31), which usually increases the current. However, this is not always the case. If the transmission through some channels which have major contribution to the current is greatly suppressed, it may not be compensated by the increase of channels. Fig. 310 clearly shows such an example, where as the bias voltage 0.0 0.2 0.4 0.6 0.8 1.0 Volts Figure 311. The differential conductance as a function of bias voltage for both trans and cis configurations. 1.0 0.8 o 0.6 E S0,4 I 0.2 0.0 6. 5 5.5 5.0 Energy (eV) 4.5 4.0 Figure 312. The transmission as a function of energy at 0.64, 0.72 and 0.8 Volts for cis configuration. ., 2 Volts \; J Q0, 2 Volts 0 Volts Ef 29 increases from 0.72 Volts to 0.8 Volts, the transmission through the HOMO state is suppressed thereby leading to the negative differential conductance in Fig. 311. CHAPTER 4 QUANTUM TUNNELING THROUGH MAGNETIC JUNCTIONS 4.1 Introduction Quantum transport in open systems far from equilibrium is an important and challenging problem [15] as mentioned in previous chapters. This problem has attracted considerable attention in the context of semiconductor based devices for which the relevant electronic structures can be modeled using the effective mass approximation. More recently, considerable progress has been made with respect to application of quantum transport theory to molecular electronics. In particular, the effect of finite bias voltage has been studied with a number of first principles approaches [620] with various degrees of success. The preceding chapter is an example. The problem of spindependent tunneling in magnetic heterostructures, however, presents a significantly greater challenge than the spinindependent case to the theory of nonequilibrium transport because one must simultaneously deal with the complexities of nonequilibrium transport and complex electronic structures along with magnetism. Spindependent tunneling junctions are considered to be the most likely candidate for nextgeneration devices for spintronics. Compared to the trilayer or multilayer giant magnetoresistance (GMR) systems, magnetic tunneling junctions tend to have a relatively higher magnetoresistance ratio [4548], which is an important factor in device designs. The magnetoresistance ratio is defined as (GpGa)/Ga, where Gp is the conductance for the case that the alignments of net magnetic moments of two leads are parallel and Ga is the conductance for the case that the alignments of net magnetic moments of two leads are antiparallel. One of the significant differences between a GMR system and a spindependent tunneling junction is that the former almost always operates within the linear response regime, while the latter usually has at least a bias voltage of a fraction of one Volt across the distance of a few nanometers, placing it far away from the equilibrium state. It has been observed experimentally [49] that the bias voltage across a tunneling junction can greatly reduce the tunneling magnetoresistance (TMR). Understanding the dependence of the TMR on the voltage bias has become one of the critical issues in the development of spin dependent tunneling devices. There have been several firstprinciples studies [5054] of tunneling conductance of trilayer junctions of the form FM/barrier/FM, where FM is a ferromagnetic material and the barrier layer is either a semiconductor or an insulator. These studies were limited to the linearresponse regime. Researchers generally found exceptionally large TMR ratios compared to experiments, and strong dependence of TMR on the interface resonance states in the minority spin channel [5152]. It is not clear whether the interface resonance states remain relevant under large voltage biases, which will significantly modify the interface electronic structure. Among the finite bias studies, Zhang and White [49] constructed a phenomenological model to explain the voltagedependence of TMR in terms of the impurity scattering at the interface. They suggested that with a high quality barrier, a different voltagedependence of TMR may occur. Some groups [5557] have proposed different models suggesting that even without impurity scattering, the TMR also decrease with the bias. One work attempted to explain the TMR reduction in terms of magnon excitations at the interfaces [58]. However, these models are often based on overly simplified scattering potentials, e. g., simple step barriers, and other assumptions such as neglecting the k// (Bloch wave vector in two dimensional plan perpendicular to the current) dependence of the transmission and reflection probabilities. Such simplifications have been shown to lead to qualitatively incorrect conclusions in linear response models [59, 60], including incorrect dependence of the conductance and magnetoresistance on the barrier thickness and height. Today, firstprinciples calculation [620] of electronic structure and tunneling conductance under a finite bias voltage has been successfully applied to molecular electronics. As mentioned in Chapter 1, the firstprinciples approach can be built on either nonequilibrium Green's functions formulation or scattering wave formulation and the scattering wave has so far been implemented only with jellium electrodes, due to the difficulties associated with computing and storing Bloch wave functions of realistic leads. The layer KorringaKohnRostoker (KKR) approach [61], which was used successfully in calculating linearresponse transport properties of spindependent tunneling junctions [51, 52], uses an extended basis set, spherical waves inside atomic spheres and plane waves in the interstitial region. The scattering wave solution for incident Bloch wave functions is provided in Ref. 51. This allowed straightforward calculation of the twoterminal conductance within the linearresponse regime. By extending this approach to treat finite bias, we can overcome the shortcoming of the scattering wave method and treat realistic leads. The test system for our calculation is a Fe/FeO/MgO/Fe(100) tunneling junction. This system has good lattice match between the Fe substrate and the MgO lattice. It has been successfully grown experimentally with demonstrated epitaxial growth [6265]. The inclusion of a single atomic layer of FeO on one of the interfaces is important. An earlier linearresponse calculation [52] for the Fe/MgO/Fe(100) tunneling junction showed a TMR ratio that is two orders of magnitude larger than experimental measurements [65]. In a subsequent calculation [66] this disagreement between theory and experiment was shown to be mostly due to the presence of the FeO layer on the interface which was observed experimentally [63, 64]. A finite voltage bias will change the interface electronic structure, which in turn will impact the TMR ratio. The Fe/FeO/MgO/Fe(100) tunneling junction would be a good system to examine whether the change in electronic structure due to the bias voltage alone can explain the change in TMR. Three effects of finite bias need to be considered in a firstprinciples theory. The first effect arises from the energy dependence of the number of channels in the wire and, in general, of the transmission probability of the conductor or the tunneling barrier. This can be accounted for by replacing the linear response formula with an integral over the energy. The second is the change in the electronic structure of the sample due to the Coulomb field of the accumulated charge induced by the flow of the current. This is usually small for good conductors, but can be significant in a tunneling device. The third one is the effect of the current on the exchangecorrelation potential in a density functional theory approach. Most forms of this potential are developed for equilibrium electron systems, and those incorporating a current are usually rather complex and difficult to implement numerically. In this work we only consider the first two effects of a finite bias, and neglect the effect of a current on the exchangecorrelation potential. In the following section, we will show that our approach is equivalent to the Keldysh Green's function approach within the singleparticle picture, but without the need for the truncation of the matrix elements in the Hamiltonian and the overlap matrix as required by other first principles methods using the Keldysh Green's function. 4.2 Nonequilibrium Green's Function Formalism in Real Space The Keldysh formalism incorporates the open boundary conditions of a non equilibrium system into the components of the Green's functions. Because not all of the Keldysh Green's functions are independent of each other, we can use the usual advanced and retarded Green's functions along with one additional component of the Keldysh Green's functions. Common practice is to use the 'lesser' Green's function G' which is defined as G' (r, t; r t) = i ( (r',t ) (r,t)) (41) where the angle brackets indicate the ensemble average and V'+, are field creation and destruction operators of electrons. There is a simple relation between the lesser Green's function and the electron density p(f, t) = iG (F, t; F, t). Since we are only interested in the steady state solution, the time dependence is only through At = t t'. So, the lesser Green's function can be Fourier transformed into G' (F, F'; At) = G' (r, '; E)e'EAtdE 1 e (42) G'(,';E) G'(F, ';t)eEtdAt Note that writing the Green's function with a single energy argument implies that the scattering in the tunneling region is elastic. The charge density is also independent of time, () = 2fG' (,,E)dE. (43) Z.7U> ' The problem remains to evaluate G' (r,'; E). For noninteracting Hamiltonians, the approach in the literature [20, 22] requires the decoupling of the two electrodes. This effectively truncates the range the Hamiltonian at a finite distance, and may not be appropriate for some systems. Here, we take a different approach using the scattering wave representation. Our approach will not be easily generalized to interacting electron systems, but is valid within the typical first principles methodology based on the density functional theory. We rewrite the Fourier transformed Green's function in the form G<(i,',E)= /' (r',t + At)V/(r,t))e'EAtdAt. (44) In the framework of meanfield theory, these two operators can be expanded as single electron operators + (ir, t)= Y0(r, t)C + n (45) y(r, = On(r, t)C, n where o, is the stationary state single particle wave function and two terms C,, C, are creation and destruction operators for this stationary state. Substituting equation (45) into equation (44), we have G< (, ',E)= i 9(E E) J Cn )CC n (r')o (r) (46) = ij 3(E En )fno (r')n (r) n where the product C C, is the number of particles that occupy the state o and its ensemble average is the distribution function f,. Equation (46) is a general expression for the lesser Green's function in terms of the stationary wave functions and the corresponding distribution function. Now, we need to adapt this expression for the tunneling junction system. An open system connected to two electrodes has a general geometry as depicted in Fig 11. As before the whole system is divided into three regions, the left lead, the device region and right lead. The wave functions and the distribution function must be considered separately within each region and continuous at boundaries. There are two types of stationary states: scattering states and bound states. The basic assumption of our approach is that the left lead maintains an equilibrium distribution function f (E /L ) which is diagonalized by scattering function '+ incident from the left lead, and likewise for the right lead distribution function f (E / ) and the scattering wave functions from the right lead 0 Both f, and f are FermiDirac distribution functions defined with their respective chemical potentials /L and /R Next, using the subscript / to denote all degenerate indices of '+ and m correspondingly foro equation (46) can be written as (neglecting the bound states for the moment) G' (r,r',E) = if, 18(E E,) 0(r' ), (r)+if 1,8(E Em)m(r')^m(r). (47) 1 m We compare this result with that obtained in Chapter 2. Rewrite equation (211) as G' = fG'FLG" +fGrFRGa. (48) Both equations are valid within the single particle picture. However, equation (47) is more general than the equation (48), as it does not require the complete decoupling of the two electrodes, thus no truncation in the range of the matrix elements of the Hamiltonian is needed. In particular, equation (47) applies to the entire space whereas equation (48) can be applied only within the device region. Substituting equation (47) into equation (43), we obtain the scattering part of the electron density p, (r) = f (E k) 2 + f (E R) 2. (49) / m Including the contribution from the bound states, the lesser Green's function contains an additional term i 3(E Eb )fb* (r) (r') where fb and b are the distribution b function and the singleparticle functions of bound states. Accordingly, the electron density also has an extra contribution from p Zf (E)b 2 (410) b The determination of fb is not unique and is still being investigated by researchers. 4.3 Selfconsistency and Tunneling Current for A TMR Junction In this section we outline our approach for nonequilibrium selfconsistent calculations of electronic structure of a TMR junction. Although the approach applies generally to any leaddevicelead ensemble, the particular Fe/FeO/MgO/Fe tunneling junction in our calculation has a layered structure with a twodimensional periodicity along each layer. This approach allows a two dimensional Fourier transform that yields a quasi onedimensional problem for each reciprocal vector k//. The left and right leads are semiinfinite in the direction perpendicular to the layers. They maintain different chemical potentials P, and /R . Without loss of generality, we may assume/ u > /R The device region includes all of the tunneling junction and a sufficient number of layers of materials from the left and right leads, so that the potential and the charge density in the lead regions can be taken to be the asymptotic values. Within each electrode, the Bloch eigenstates can be divided into rightgoing states (+ and leftgoing states ( The scattering states when there is a device region can be categorized into those incident from the left lead traveling to the right, !+, and those incident from the right lead traveling to the left, 0! Their wave functions can be expressed in terms of the Bloch eigenstates within each electrode, and a scattering wave function in the scattering region, in the form p + tl ,l, Left lead = (+, device region (411) to (0' Right lead m and ( +t ,t' ,,,(+', Right lead S= ( device region (412) t'lm( Left lead where t, t' are the transmission coefficients and r, r' are the reflection coefficients. The superscripts + indicate the direction of travel along the z direction, and the indices 1, m label Bloch states in the left and right electrodes, respectively, and span over all possible values of k// and the spin. If one takes the jellium approximation in both leads, the Bloch wave functions become plane waves and 0+ become the scattering states used in Ref. 10. With the help of the layer KKR approach, we have gone beyond the jellium approximation and are able to deal with the realistic leads. At zero temperature, both distribution functions become the step functions (E /L, ). Then from equations (49) and (410) we obtain p(r) = dEp (r) + dEp (r) +pb(r) (413) where p =(r) YZ (E E,) , 2 (414) p (r) = Z (E Em) ( m These equations hold for all regions including the device region and both leads. We can rewrite equation (413) as p(r) = p (r) + dEp (r) (415) PlR where the equilibrium part of the electron density corresponding to '/R is defined as PlR p (r)= dE[p+(r) +p(r) + pb(r) (416) Note that the superscript R indicates the chemical potential /RU used in the calculation of Pe, yet the density pe is defined for all three regions. Now, let us consider the electron density in the two leads separately. In the right lead, the first term in equation (415) should give the asymptotic value of the electron density far away from the device region at zero bias. The second term, the non equilibrium part, is due to the electrons transmitted from left. Neglecting the interference terms between different Bloch states, which should be averaged out in an ensemble average, we have PE within the right lead p(r)= 3(E Em)tm (r) (417) Im The electron density in the left lead contains more terms than the right lead. Again, neglecting the interference terms, we can write p (r)= c(E E,) p (r)2 + 3(E E,)R, R1,, (r) (418) 1 11' The normalization of the Bloch wave functions is such that the integral over the volume of a unit cell d3r L (E E,) 9, + 5(E E,) p equals to the density of states at energy E of the left lead, and likewise for the right lead. The selfconsistent procedure for the nonequilibrium calculation is very similar to the usual equilibrium calculation for an interface system, except in this case, the charge densities and the electrostatic potentials for all three regions need to be iterated simultaneously, unlike the equilibrium case in which the semiinfinite (lead) parts are calculated using the bulk selfconsistent calculations before the interface region is calculated. An additional consideration is given to the relative shifts in the electrostatic potential between the three regions. Although the leads should screen any net charge quickly so that the overall system should be charge neutral, the effect of dipole layers which cause constant shifts in the potential cannot be screened. The size of the dipole moment on each layer also depends sensitively on the convergence of the calculation. To overcome this problem, separate constant shifts in the potential are added to the Kohn Sham potentials for the three regions. These shifts are adjusted at each LSDA iteration until the net charge within each region reaches zero respectively. In general the shifts in the potential differ from the bias voltage, as has been shown in earlier studies [21]. In order to see how the spindependent calculation is done, we write the electron density in equation (414) in terms of spin explicitly PE,(r)= Z8(E Ek) 2 (419) P.,(r) =3(E Ek, O) k, where the term a is the spin index. The subscripts l,m now comprise both k// and a . Remember that PE+ sum over ,1 e a and p, sum over p,,m e u. Then equation (4 15) can be written as p,(r)= p", + fdEpE (r) (420) PlR where the equilibrium part is pq, (r) = JdE[p (r) + PE (r)]+ p, (r). (421) The last term means only those bound states which have spin c are counted. The equation (421) can be used to calculate the electron density for each spin over all space. After obtaining the electron density for each spin, we use the common local spindensity functional to calculate the effective potential. Similarly, the current density for each spin is given by a Landauertype formula J = eJdE e ti m (422) m,lJc Vl where v,, v, are the group velocities of the Bloch states, and only the Bloch states for the a spin are summed. The velocity factors are necessary because the transmission coefficients are defined based on the continuity of the wave function, not the flux. The equilibrium term in equation (420) can be evaluated using the standard contour integration technique for Green's function methods. The nonequilibrium term on the other hand, is more difficult to evaluate numerically. The integration over the energy for the nonequilibrium term has to be computed along the real axis, which is numerically noisy due to the singularities and resonances arising from the threedimensional Bloch states in the leads. This problem is not unique to the scattering wave function approach employed here. Even if one were to compute the lesser Green's function directly, the integration still would need to be done along the real axis because the lesser Green's function is not analytic in both upper and lower complex energy planes. Thus a complex contour cannot be used in either approach to reduce the numerical noise of the non equilibrium term. The numerics of the real axis integration can be improved greatly using the technique by Brandbyge et al [33]. We note that equation (420) can be equivalently written as I/' ', (r)= p dEp, (r) (423) PlR where /PL P~(r)= dE [p,(r)+p p(r)]+,p (r). (424) In principle, the electron density evaluated from either equation (420) or equation (424) should be identical. But because these equations use two different energy contours, the numerical errors are different. It is shown that the numerical error is minimized by combining the two electron densities in the following manner [67]: p'W = CP + (1 c) p' (425) where L 2 c neq,T (426) SL R [neq, +[pneq,a] with ,L = f dEp Pneq,T d ,( PR (427) = f dEpE, PR 4.4 Application to the Fe/FeO/MgO/Fe Tunneling Junction The Fe/FeO/MgO/Fe junction can be grown by depositing MgO epitaxially onto Fe whiskers and then depositing another Fe electrode epitaxially onto the MgO layer. Experimentally it was shown that a single atomic layer of FeO forms when MgO is deposited onto Fe [6264]. For the MgO/Fe interface, we use the same structure as in Ref 52. The lattice constant of the Fe layer is fixed at 2.866 A. In order to match [100] layers of two materials, the MgO lattice constant is taken to be a factor of /2 larger than that of Fe. The distance between Fe and O layers is taken to be 2.16 A, while the distance between Fe and O atoms in the FeO layer is taken to be 2.1 A. Because we are unable to incorporate the considerations of chemical or structural disorder into the Keldysh Green's function formalism, a necessary step to account for the under occupancy of the oxygen site in the FeO layer, we assumed 100% oxygen in our calculation. Numerical convergence in terms of the number of energy points between two Fermi energies and the number of k// points in the twodimensional reciprocal space needs to be considered carefully because the integration is along the real energy axis, and can be prone to be numerical errors, as just discussed. For the selfconsistent part of the calculation, we use one energy point for biases below 0.01 Hartree and two energy points for biases between 0.01 and 0.018 Hartree, the largest bias in our calculation. We find increasing the number of energy points further has little effect on the selfconsistent solution. For tunneling conductance calculation, we begin with three energy points in a uniform mesh for a bias of 0.001 Hartree. Then for each increment of 0.001 Hartree in the bias voltage, two additional energy points are added to the mesh. Our calculation at the bias of 0.001 Hartree showed that, when the number of k// points is increased from 2800 to 8262, the change in the conductance of parallel alignment of the magnetic moments in the electrodes is less than 1% and the change for antiparallel alignment is less than 3%. Therefore, 8262 k// points are used in the integration for all of the conductance calculations. In the following discussion, the majority and minority spin channel always refer to the left Fe lead. 0.3 Majority + Minority 0.2 0 QMg > 0.1 (1D a, 0 S 0 Fe FeO Fe c 0.1 MgO + 0.2 0.3 Figure 41. Shifts in the density of states in each atomic layer at a bias voltage of 0.54 V. It is interesting to examine the change in the electronic structure due to the applied voltage bias. The main effect of the voltage bias here is the accumulation of the charge on the interface layers next to the barrier. We find that the change in the net charge on the interface layers is small, and this change does not cause significant changes in the profile of the electronic density of states at each layer. We see the nearly constant shifts in the energy in the barrier that should be consistent with the shift in the electrostatic potential at each layer. The shifts in the DOS between E=0.5 and E=0.5 Hartree for each layer for the bias voltage of 0.54 V is shown in Fig. 41 for parallel alignment of the moments. The DOS for the iron layers in the leads is shifted by roughly +1 / 2 V on either side of the barrier, and the shifts for the barrier layers correlate monotonically with the layer positions. The shifts for both spin channels are about the same for each layer. Two interesting features occur near the interfaces. On both interfaces, the first insulator layers (FeO or MgO) have the same shifts as the electrodes. This is somewhat surprising, noting that such alignment does not conform to the extrapolation of the shifts of other MgO layer. The nearly exact alignment of the DOS of the interface insulator layers with the electrodes does not agree with the model given by Ref 34 in which an oscillatory behavior is predicted for the electrostatic potential near the interfaces. On the other hand, we find that the shifts in the DOS scale essentially linearly with the bias voltage, as one would expect. For negative biases, the shifts are the same in magnitude with the opposite signs. From the linear scaling of the DOS shifts with the bias voltage one would expect that the capacitive charge should also scale linearly with the voltage. Surprisingly, this is not the case. In Fig. 42 we show the capacitance dQ / dV as a function of the voltage. There are three different charges plotted in this figure, none of which strictly conforms to the usual definition of capacitive charge, due to the difficulty of separating the capacitive 0.01 Left electrode 0 0.008  O o 0.006 Right electrode , >x x," S0.004  S,+ Net charge C 0.002 0 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 Bias Voltage (volts) Figure 42. Derivative of the capacitive charge with respect to the bias voltage. The charge on the right electrode includes that on the oxygen atom in the interface MgO layer. 'Net charge' includes that on both the electrode and the dielectric charges on the MgO layers on the same side of the midpoint of the barrier. charge in the electrodes and the dielectric charge in MgO. The capacitive charge must be antisymmetric, i.e., with Q on opposite sides of the barrier. This antisymmetry is not possible with any partition of the atoms in our calculation, unless we include the charges on all iron and MgO layers on the same side of the midpoint of the barrier. The curve labeled 'net charge' is this sum per twodimensional (2D) unit cell. This sum is the same (with opposite signs) for both sides of the barrier. However, the sum includes both the capacitive and the dielectric charges. To separate out the capacitive charge at least approximately, for the left electrode we sum the charge for all iron layers but excluded the FeO layer. For the right electrode we sum the charge for all iron layers plus the charge of the oxygen atom on the interfacial MgO layer. The reason for including the oxygen atom is that it represents more than a quarter of the net change in the charge on the right side, consistent with our observation that this oxygen atom seems to be 'short circuited' with the iron electrode. The two curves in Fig. 42 labeled 'left electrode' and 'right electrode' represent this approximation of capacitance. For comparison a model of capacitance per unit cell for two plates separated by a vacuum gap of 21 A (which is the distance between two iron electrodes), is about 0.0044 electrons/V. Thus our self consistent calculation yields an effective dielectric constant for the MgO film (coated with an atomic layer of FeO on one side) of around 2. Since in our calculation no lattice relaxation is included, this number should be compared to the optical dielectric constant measured experimently. The value derived from the index of refraction measurement [35] (n=1.7) for bulk MgO is s = n2 =2.9. 0.04 . 0.035 Right 0.03  o 0.025 0.02 Left E 0.015  0.01 Bulk electrode (left) 0.005   0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 Bias Voltage (volts) Figure 43. Derivative of the moment with respect to the bias voltage. The labels 'left' and 'right' indicate that the layers to the left (right) of the midpoint of the barrier are summed. The contribution to the sum from the nonbulk layers of the iron electrodes is taken to be the difference between the moment of that layer and the bulk moment in the electrode on the same side. Because of the neutrality of the total system, the total net charge integrated over the 2D unit cells of all the layers on one side of the midpoint of the barrier should be the same with the opposite signs from both sides. This charge is shown as 'net charge' in Fig. 42. Evidently, it is not symmetric with respect to the bias voltage. For positive biases, it drops quickly as the voltage increases. At about 0.5 V, it becomes negative. This means that at this voltage, any increase in the charge in the electrode due to the increasing bias is completely matched and canceled by the opposite dielectric charge from the MgO layer. An unexpected result is the change in the magnetic moment at the interfaces for parallel alignment of the moments. For a symmetric junction such as Fe/MgO/Fe, if the net moment increases with bias on one interface, it must decrease on the other interface due to symmetry. For the asymmetric junction ofFe/FeO/MgO/Fe, we find that the moments on both interfaces increase with positive bias and decrease with negative bias. Thus the derivative of the net moment with respect to the voltage is always positive, as shown in Fig. 43. Unlike the charge, there is no indication of saturation of the moment changes at the interface as a function of the bias. The asymmetric changes in the moments mean that the majority spin electrons react differently to the bias voltage than the minority spin electrons. Since the rate of the change in the moment is almost an order of magnitude larger than the rate of change in the charge, charge accumulation in the two spin channels must be opposite in sign. The curve labeled 'bulk electrode' in Fig. 43 is the change in the moment in the bulk for the left electrode. The changes in the moment in the bulk regions on the two electrodes are the same in magnitude but opposite in sign. The change in the interface magnetic moment in response to the bias voltage is detectable experimentally if one can keep the right electrode sufficiently thin. The right electrode (which does not have an FeO layer at the interface) is the last deposited iron layer and its moment changes in the opposite direction as the interfacial moment as a function of the bias voltage. At zero bias, the conductance is given by G, =e2 2V m 2 (428) Imeo V It is important that our calculation from the finite bias approach extrapolates to the linear response results obtained in earlier works [51, 52, 66]. The k//resolved conductance at zero bias, shown in Fig. 44, reproduces all of the features seen in earlier works. For the parallel case, the majority spin channel is dominated by the contribution from a region around k//=0. The conductance due to the minority spin channel is about two orders of magnitude smaller and is not shown in the figure. For the antiparallel case, our calculations show that the predominant contribution to the minority spin current comes from tunneling through the interface resonant states which correspond to sharp peaks of the transmission. The calculation gives a negative TMR of 28% at zero bias. We note that an earlier work obtained a small positive TMR ratio for the same junction [66], but our calculation yields a small negative TMR ratio. This difference appears to be due to the extreme sensitivity of the minority spin channel conductance to the convergence of the integration over k//. The earlier calculation used 2800 k/ points while we used significantly more k/ points (8262). Nonetheless, the general trend studied, that of TMR as a function of oxygen concentration in Ref. 31, and the TMR as a function of bias voltage studied here, should remain valid as long as the number of k// points used is consistent within each study. Transmission (105) .0 0.2 .5 0.15 S0.05 00' Transmission (104) 4.5  4.0 I 3.5 I 3.0 I 2.5  2.0 J 1.5  1.0  'II I I. 0.5 " .. "'^. .: .. ..Z  0.2 0 15 0 05 Ky Figure 44. Zero bias transmission probability as a function of k//, (a) Majority spin channel for parallel alignment of the moments; (b) The dominant spin channel for antiparallel alignment of the moment. 40 L 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 Bias Voltage (Volts) Figure 45. Current density as a function of bias voltage for majority and minority spin channels for the case of parallel alignment of the moments between electrodes. 6 4 .'Cu 2 0 2 4 6 8 . 10 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 Bias Voltage (Volts) 0.3 0.4 0.5 Figure 46. Current density as a function of bias voltage for majority and minority spin channels for the case of antiparallel alignment of the moments between electrodes. The tunneling current is calculated through equation (422). The effect of the bias voltage, through the change of the electrochemical potentials in the two electrodes, enters this equation in two ways. First, the bias voltage can change the transmission probability for each conduction channel, i.e., for each k// point at each energy. Second, when the bias increases, the number of conduction channels in the leads that contribute to tunneling also increases. An interesting question is how much of these effects is captured by rigidly shifting the electrostatic potential for each atomic layer and neglecting the effects of charge rearrangement, as was commonly done in many early nonself consistent calculations, and, conversely, how important is charge selfconsistency in accurate determination of the TMR. For this purpose, we compare the fully self consistent calculations of the tunneling current and the TMR with a calculation constructed from a selfconsistent zero bias potential, in which the electrostatic potentials in the left and right electrodes are rigidly shifted by + 1/ 2 V where V is the bias voltage, and the potentials of the barrier layers are shifted linearly with the bias according to their positions. Because of the asymmetry of the tunneling junction, the positive and negative biases give different tunneling currents. Figure 45 shows the IV curves for the majority and minority spin channels for parallel alignment of the moments. As expected, we see a large linear regime in this figure. For low biases, the slope of the curve gives the linear response conductance a value ofl. 15 x 107 /(m 2). This number agrees very well with the result obtained from equation (428) which is1.17 x 107 /(Qm2) Note also that the nonself consistent calculation yields the same results within the linear region (for bias voltage up to about 0.1 V) but deviates significantly at higher biases, reaching about twice the selfconsistent results at 0.5 V. The current of minority spin channels actually 3 .5 . '0.0_Volts' '0.054 Volts' '0.108 Volts' .......... S2.5 0 c 2.0 0 E 1.5  1.0 0.5 0.2 0.15 0.1 0.05 0 0.05 0.1 0.15 0.2 Kx (Ky) Figure 47. Transmission probability at energy E = /R as a function of k/ along the line of kx = ky for the majority spin channel for parallel alignment of the moments at different bias voltages. shows an oscillatory behavior. Since it is negligibly small compared to the majority spin channels, we cannot see this in Fig. 45. Figure 46 shows the total tunneling current for antiparallel alignment of the moments. The curve labeled 'SCF' is in fact calculated by simply flipping the moment alignments of half the layers of the corresponding parallel SCF calculations, thus this is not a truly SCF calculation for the antiparallel alignment. However, one fully SCF calculation for the antiparallel alignment was performed at V=0.5 V and the difference is negligible both in terms of charge transfer and tunneling current. The curve labeled 'non SCF' is calculated with the rigid shift of a selfconsistent zerobias electrostatic potential. For most of the bias voltages the tunneling currents agree between the 'SCF' and 'non SCF' calculations, except for large negative biases. The tunneling current exhibits an oscillatory behavior as a function of positive biases, probably due to the matching and mismatching of interface resonance states on both sides of the barrier as the bias voltage moves these states in energy and k/. We will see this later from the analysis of integrated transmission. Let us consider a positive bias as an example to study the effects of bias on the tunneling process. By positive bias we mean that the electron flux is from the left (FeO side) to the right, or that the current is from the right to the left. We calculated the transmission probability for the spin channel with the largest contribution to the current at the energy E = /u and along a line of kx = k, in the 2D reciprocal space, under different positive bias voltages. For the case of parallel alignment of the moments in both leads, Fig. 47 shows the results of the majority spin channel at bias voltages of 0.0, 0.054, and 0.108 V. All curves exhibit very similar structures. It is clear that the bias does not affect the transmission very much. This suggests that the linear regime of IV curve for the majority spin channel may extend up to 0.108 V. For the antiparallel case, Fig 48 shows the results for the same range of bias voltages. Although all curves has the similar overall structure (the transmission is low at k/=0 and have one peak on each side of k/=0), the position, height and width of these peaks change dramatically as functions of bias. From 0.0 to 0.082 V, the peaks are enhanced, and from 0.082 to 0.108 V, the peaks are reduced. As discussed in an earlier linear response study [20], the position and the strength of these peaks are determined by several factors. Among them are the amplitude of the wave function at the interfaces on both sides of the barrier, and the decay rate across the barrier as determined by the complex band structure in the barrier. When interface resonance states exist, the amplitudes of the wave function can be enhanced by several orders of magnitude at the interface through coupling to those interface resonance states. This coupling gives rise to the sharp peaks in the transmission probability at certain k// points. Applying a bias voltage or changing it will change the energy and k/ positions of these interface resonance states, in particular the relative positions of interface resonance states on both sides of the barrier can be moved. This in turn can cause large changes in the transmission probability at corresponding k// values. 3.0 , '0.0 Volts' '0.054 Volts'  2.5 '0.082 Volts' '0.108_Volts' i o 2.0 \ 1.5  E i c 1.0 C0 0.5 ; 0.2 0.15 0.1 0.05 0 0.05 0.1 0.15 0.2 Kx (Ky) Figure 48. Transmission probability at energy E = /R as a function of k/ along the line of kx = ky for the minority spin channel for antiparallel alignment of the moments at different bias voltages. As mentioned at the beginning of this section, when the bias voltage increase, the number of conduction channels in the leads that contribute to tunneling also increase. In order to examine this effect, we define the integrated transmission as 1 " T(k,) t(k, )dE. (429) S R /'R In Fig. 49 and 410, we show this integrated transmission probability as a function of k// along a line of kx=ky in the 2D reciprocal space. The patterns shown in these two figures closely resemble those in Fig 47 and 48. However, the integration over the energy has smeared out to some degree the effects due to the interface resonance states, so that the change in transmission probability as function of bias is less dramatic. Based on these figures, we can expect that the majority spin tunneling current for the parallel case should scale approximately linearly with the bias, and the minority spin tunneling current and the antiparallel tunneling current should increase more slowly than linear with the bias. This behavior should lead to an increase of the TMR ratio as a function of the bias voltage. 3 .5 . '0.0_Volts' '0.054 Volts'  U 3.0  3.0 '0.108 Volts' .......... c 2.5 C,/) .(n 2.0 C S 1.5 0 S 1.0  0) 0.5  0.2 0.15 0.1 0.05 0 0.05 0.1 0.15 0.2 Kx (Ky) Figure 49. Integrated transmission probability as a function of k// along the line of kx = k, for the majority spin channel for parallel alignment of the moments at different bias voltages. Although the negative bias affects the system differently, the general feature is the same: the transmission probability of dominant channel of antiparalle case is much more sensitive to bias than that of parallel case. At small bias, the antiparallel tunneling current also increases more slowly than the parallel tunneling current as a function of bias voltage. Figures (412) and (413) show the integrated transmission as a function of k// for all k// points at a higher bias of 0.378 V. Compared to Fig 44, the broad peak in the parallel case becomes broader, a result of the integration over the contribution from all channels at different energies. For the antiparallel case, the structure of the peaks due to the interface resonance states becomes more complicated than at zero bias. The TMR at finite bias is defined as (J J) / Jp, where Jp, J, are total currents for the parallel and antiparallel case respectively and its linearresponse limit is (Gp Ga)/Ga as given in the introduction. The voltage dependence of the TMR is shown in Fig. 413. Previous calculations that used simplified models [49, 50, 58] yielded TMR ratios that are large and positive at zero bias, then decrease rapidly with bias. Our calculations yield a very low TMR ratio at low biases, even negative at 0 and + 0.027 V. Then, for positive biases, the TMR ratio increases rapidly as the bias voltage is changed from 0.054 to 0.108 V. This increase is caused by the decrease of the antiparallel current. When the bias voltage is larger than 0.108 V, the TMR ratio starts to oscillate with the bias, but overall still rises with the bias, albeit more slowly. For negative biases, the TMR increases from 0.027 to 0.432 V and decreases at 0.486 V, which is caused by slow increase of parallel tunneling current. The comparison between the SCF calculation and the nonSCF calculation for TMR is similar to that of the tunneling current for parallel spin alignment. This indicates that the charge selfconsistency is more important for the majority spin channel than the 0 1 1I 1, 0.2 0.15 0.1 0.05 0 0.05 Kx (Ky) (a) 0  0.2 0.15 0.1 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2 Kx (Ky) (b) Figure 410. Integrated transmission probability as a function of k// along the line of kx = ky for the majority spin channel for parallel alignment of the moments at different bias voltages. Transmission (105) 0.5 Kx 0.1 0. .0I S0.2 0.15 ~ 0.1 ' 005 Ky Figure 411. Integrated transmission probability as a function of k// for the majority spin channel for parallel alignment of the moments at bias 0.378 V. Transmission (105) 1.4  1.2  1.0 I 0.8  I.1 0.6 . 0.4 0. 0.2 15 00. i *C ..i : z .o Ky Kx Ci^i i I 0 2 1 Figure 412. Integrated transmission probability as a function of k// for the minority spin channel for antiparallel alignment of the moments at bias 0.378 V. minority spin channel, presumably due to the sensitivity of the majority spin tunneling current to the chemical bond configuration at the interfaces [66]. The tunneling through the interface resonance states leads to the sharpest peaks of the transmission. Due to the presence of roughness and impurities at the interfaces, the interface resonance states may not be present or sharply located in experimental measurements. In Fig. 413 we also indicated the TMR ratio at zero bias if we removed all of the contributions from the interface resonance states. Although the TMR ratio increases significantly, this increase does not change the overall trend of increasing TMR ratio with bias. Thus we conclude that at least for Fe(100) electrodes, where a single atomic layer of FeO is formed at one of the interfaces, the TMT ratio is close to zero at small biases, and should increase with the bias voltage. We believe that this conclusion should also apply to barrier materials other than MgO, especially when the barrier layer is asymmetric as a result of the order of deposition, so that the TMR ratio is low at small biases. This conclusion supports the suggestion by Zhang and White [9], that for high quality barrier, the voltagedependence of TMR may be different than the observation in earlier experiments that the TMR ratio decreases with bias. For the system we studied, this unusual behavior is based on following facts. First, for the parallel case, the transmission of dominant conduction channels (peak for majority spin channel) is not sensitive to the bias. This gives a large linear regime as a function of bias. For the antiparallel case, the dominant conduction channels are sensitive to the interface electronic structure, in particular, interface resonance states, which in turn are sensitive to the bias voltage. Often the maximum of transmission occurs when a tunneling state couples efficiently to both leads. This may be the case when interface resonance 'TMR:SCF' B i 3.5 'TMR:NonSCF'  3 2.5 1.5 \ 1 . 0.5 0 0.5 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 Bias Voltage (Volts) Figure 413. Magnetoresistance ratio as a function of bias voltage. Point A denotes the MR ratio at zero bias by removing all contribution from interface resonance states. states are present on both interfaces. When the bias increases, it is likely that the once efficient coupling to both leads is greatly reduced due to the change in the interface resonance states. This would lead to a great reduction in the tunneling current for the minority spin channel and for the antiparallel case. Thus, it is possible for TMR to increase with the bias. Although one may expect that for amorphous barriers the effects of interface resonance states are much reduced, the impurity states at the interface may play the same role as the interface resonance states in the system we discussed in this work, so that the argument presented here may still be qualitatively true. Thus, it may be possible to observe the similar behavior of TMR even for amorphous barriers. 4.5 Conclusion In conclusion, we have completed the first spinpolarized calculation using the nonequilibrium firstprinciple method for spindependent tunneling junction under a finite bias. Application of this method to the layer KKR treatment of the Fe/FeO/MgO/Fe tunneling junction leads to surprising results that are not anticipated by simple models but are reasonable in view of the change in the electronic structure at the interface in response to the bias voltage. Our calculations indicate a negative TMR ratio at very low biases and a rapid increase of the TMR ratio at finite biases. We believe that this behavior is not limited to the Fe/FeO/MgO/Fe tunneling junctions and may be observed experimentally in other systems. This result has import practical implications since increasing the TMR ration at finite bias voltages is one of the important goals of spintronics research. Our calculations also show significant deviations from nonSCF calculations using rigidly shifted electrostatic potentials at bias voltages about 0.1 V, as well as the nonlinear capacitive charging effect and the asymmetric magnetic moment change as a function of bias voltage. All these results highlight the important role of a fully selfconsistent nonequilibrium calculation. CHAPTER 5 PROBLEMS IN THE CURRENT NONEQUILIBRIUM APPROACH FOR OPEN SYSTEMS BASED ON DFT 5.1 Introduction As discussed in previous chapters, for the open system under a finite bias, the non equilibrium firstprinciples approach has been proved to be successful. The approach is based on density functional theory, and the typical procedure can be summarized in two steps as we mentioned in Chapter 2: the first is the selfconsistent procedure based on DFT to determine the effective singleelectron potential. Here, of course, the density functional used is valid for ground states only, but is assumed to be valid for non equilibrium systems because no alternative approaches exist. The second step is to use the Landauer conductance formula to calculate the current. The Landauer formula relates the current to the integral of the singleparticle transmission function over energy. Such approaches have been used in many systems with varying degrees of success: the shape of the IV curve can be compared to the experiment and the structure of IV curve can be explained by molecular orbitals. All these achievements are important for the understanding and future design of molecular electronics. However, large discrepancies are often seen between theoretical calculations and actual measurements of IV curves. Among possible sources of error in the theory, the use of the particular type of the density functional in a calculation is often a cause of concern. Much progress has been made in finding an energy functional form that works well for both the metallic leads and the molecular junction. Unfortunately, the problem is more fundamental than just the functional. The basis of DFT is the first HohenbergKohn theorem [23], which states that for the ground state, the external potential Vex,(r) is determined, within an additive constant Vo, by the electron density p(r). Since all ground state properties are determined by the Hamiltonian, which is uniquely defined once the external potential is determined, the onetoone correspondence between Ve, (r) and p(r) guarantees that all ground state properties are determined by the electron density p(r) In other words, p(r) is the basic variable of the system at ground state and all physical quantities can be expressed as functionals of p(r) The constant potential Vy is trivial for ground state, since if Vy changes, the Fermi energy and the electronic structure of the system simply shifts accordingly. The theorem can be easily proved through the fact that the ground state has the lowest energy. However, the same proof does not apply to an open system under a bias voltage, since the system is in a nonequilibrium state, and there is no minimum energy principle for such a state. Before we use the DFT in such a system, we have to question whether the electron density p(r) remains a fundamental variable that uniquely describes the system. In the following section, we show that the electron density alone is not sufficient to describe a nonequilibrium open system. Then we propose a modification of the HohenbergKohn Theorem for an open system at finite bias. We follow that with a numerical study that gives a functional form of the exchange energy for nonequilibrium open systems. 5.2 HohenbergKohn Theorem for Open System under Finite Bias To make our discussion more concrete, we choose a model open system with the general structure shown in Fig. 11. The system is divided into three regions: the left lead, the scattering region (the device region), and the right lead. The scattering region consists of the left contact, the barrier and the right contact. The two contacts are parts of leads that connect to the middle barrier. The leads are connected with two reservoirs at infinity (not shown in Fig. 11) and the chemical potentials of the two reservoirs are constants. A discussion of the contact between leads and reservoirs can be found in reference 68. The two leads are often to be assumed to be internally at thermodynamic equilibrium [20], but not in equilibrium with the contacts or the scattering region through the contacts. The electronic structure of the two leads far away from the scattering region are assumed to be identical to their bulk values [20], and two separate calculations based on firstprinciples method are often performed to obtain properties of the leads. Furthermore, the temperature of both leads is assumed to be zero. The bias voltage is defined as the difference between Fermi energies '/, and '/R maintained by the two leads connected to two reservoirs, respectively. The purpose of the study is to calculate the current through the system as a function of the bias voltage. Let us now examine the question whether the electron density p(r) can uniquely determine the transport properties of the system, specifically, whether p(r) alone is sufficient in determining the current through the system. Because the current is a function of the bias voltage, if p(r) is the basic variable of the system in the same way it is for a ground state, knowing p(r) should allow one to determine the voltage bias across the scattering region. In both leads, the potential and the electron density can be approximately taken as the ground state bulk value, due to the screening by the two contacts. This is the approximation generally used in the previous studies. Suppose the first HohenbergKohn theorem for the ground state can apply at least in the two leads. Then, according to the theorem, the Fermi energy /Lu is decided by the electron density in the left lead PL (r) to an constant VL and the Fermi energy /Ru is decided by the electron density in the right lead PR (r) to another constant VR. Since both leads are connected to a battery, these two constants, VL and VR, should not be independent. However, according to the ground state HohenbergKohn theorem, if we shift Fermi energies of left and right leads, thus changing the voltage bias across the scattering region, the electron density in the two leads should remain unchanged. In another words, these two constants, VL and VR, are independent if only given electron density. The result is that the voltage bias /L /i is completely arbitrary, if one is only given the electron density in the two leads. Since the bias voltage should not depend on the scattering region, we come to the conclusion that for the open system under a finite bias, the electron density p(r) alone is not sufficient to describe all properties of the system, thus conventional DFT is not applicable. The HohenbergKohn theorem must be modified in order to be applied to such systems. 5.2.1 Voltagedependence of the Density Functional In order to see how the HohenbergKohn theorem should be changed, we start from a rigorous approach to this problem based on the adiabatic timedependent procedure. We take the same timedependent process as that in reference 24 and reference 25: at time t = oo, the unperturbed system consists of three uncoupled regions, the left lead, the scattering region, and the right lead. Each one maintains its own equilibrium state (ground state for zero temperature). The Fermi energies for the left and right leads are /u and /Ru respectively. The Fermi energy of the scattering region does not need to be the same as left or right leads, but for simplicity, we assume that it is in the same equilibrium state as the right lead. The bias voltage is defined as/u Ru The couplings between leads and the scattering region turn on slowly from t = co tot = 0, such that at any time, the system is in the stationary state. The goal is to relate the current at t = 0 with the couplings completely turned on to the unperturbed Fermi energies '/L and /R Given the initial state, the properties of the timedependent system can be obtained by solving the timedependent Schrodinger equation .0 S (t) = H(t)(t). (51) The Hamiltonian can be written as H(t) = H +V(t) IH,= 1V,2+i + Vext(r,)+Z Vee(rj) (52) 2 1 I where, in the Hamiltonian, the first term H, consists of the kinetic energy, the nuclear electron interaction Vex and the electronelectron interaction V,. The second term is an extra timedependent external potential, where V can be chosen to be two barriers in the two contacts and zero in all other regions. The two barriers should be high enough such that at t = oc, the transmission probability of electrons tunneling through these two barriers is negligible; thus, the two leads and the middle barrier are well separated. The constant c should be positive and small enough such that the couplings between leads and the middle barrier turn on adiabatically. For such a timedependent system, Runge and Gross [69] showed that the wave function at time t, ry(t), can be determined by the timedependent electron density p(t) and the initial state, up to a timedependent phase factor, ,(t) = e 'O't)A/P(O, /(o0)] (53) where, [4p, r/(co] is a functional of the electron density and the initial state. Particularly, if the initial state is the stationary ground state at t = 0o, then according to the first HohenbergKohn theorem for ground state, it is completely determined by the initial electron density; hence ((oc) can be eliminated in the equation (53), so that y/(t) is determined up to a phase factor by p only. Neglecting the phase factor, we come to the timedependent DFT in this case, which claims that any timedependent physical quantity can be expressed as the functional of timedependent electron density. However, for our case, at t = oc, the system is not in the ground state. The initial state V,(oc) should be constructed to represent three local equilibrium states for three uncoupled regions with two different Fermi energies /L and /R which actually is an excited state at zero temperature characterized by /L and /R Apparently, the initial state is determined by the Hamiltonian and these two Fermi energies. Since a trivial constant added to the Hamiltonian will not affect the initial state, but will shift /'L and PR by the same amount, we know that the initial state is determined by the Hamiltonian and the voltage biasVb = /u R . Applying the first HohenbergKohn theorem for the ground state to three regions which are in local equilibrium states, the Hamiltonian can be determined by the initial electron density to a constant. As we discussed before, the bias voltage is independent of the electron density, so the initial state is a functional of electron density p and the voltage bias Vb. Att = 0, the Equation (53) can be written as, f(t) = e p(t)/[P(t),1 ] (54) which means that at a given Vb, ry(t) is determined by p to a timedependent phase factor. y/lp, b ] can be equivalently written as i/[Pb ] \ which means that the density functional form io] depends on the bias voltage Vb. We come to a conclusion that all steadystate properties at t = 0 can be expressed as functionals of p at a given bias voltage V The most important of these functionals is the exchangecorrelation energy functional Ex, []. The voltagedependence of those density functionals actually comes from the voltagedependence of the initial state. The first HohenbergKohn theorem for the equilibrium state is just a special case forVb = 0. The voltagedependence of the density functionals introduce substantial difficulty to the implementation of the theory, since we need different functional forms for different bias voltages. It is even worse than that: since the voltagedependence of the initial state varies with the system, the voltagedependence of the density functionals varies with the system, which means we don't have universal density functional for all systems even when they are under the same bias voltage. This disadvantage is caused by the fact that the basic variable of the functional, the function p(r), can not determine the bias voltage Vb, and is not sufficient to describe the system. In order to get rid of this disadvantage, we need to introduce other functions as variables of the functional, which together with p can completely determine the system. 5.2.2 The Nonequilibrium Electron Density The key to eliminate the voltagedependence of the functional in the equation (54) is to introduce other functions into the functional which can determine the bias voltage Vb and thus determine the initial state y,(oo) together with the electron density . For the nonequilibrium system, the electron density can be expressed as [20, 22], p(r)= I'd (r, E) (55) 21i where, G'(r, ) is the nonequilibrium lesser Green's function defined for the stationary state. For the system under study, which is characterized by two Fermi energies /, and /Ru (suppose /u > /,R), the electron density can be divided into the equilibrium part peq and the nonequilibrium part p,, p(r)) (r)+ pn(r) (56) where, 1 ^R 2ii P,, (r)= 2IP dsG(r,E) 7 (57) p (r) = I dsG (r, E) The term Pe, is called the equilibrium density in the sense that it can be equivalently written as an equilibrium integral [20], PLR pe (r)= Im dsG (r, E) (58) co where, GR is the usual retarded Green's function. We want to show that if given both Peg and p, at the initial time, the bias voltage is determined and then the initial state is determined. When t = oo, the left lead is well separated with respect to other regions. The left lead is in a local equilibrium state with the Fermi energy PuL while the scattering region and the right lead are in a local equilibrium state with the Fermi energy pR In this case, pn in the right lead is zero and in the left lead it also can be written as an equilibrium integral, p,(r) = Im d~GR(r, ). (59) i" PLR Since the left lead is in local equilibrium, the first HohenbergKohn theorem for equilibrium state applies. We know that in the left lead, the Hamiltonian is determined by p(r) = pe (r) + p, (r) to an arbitrary constant. If a constant V, is added to the Hamiltonian, GR (E) changes to GR (E + V0) Then, with p,, given, from equation (58), we see that PRt changes to PR + Vy. Introduce the new pR and GR (E) into the equation (59), we see that with the given p,, PL becomes /L + V0. Hence the bias voltage Vb = UL /UR is completely determined without any arbitrary constant. Then, the Equation (54) can be rewritten as, Y(t =0) = e '(t =o[pq (t = 0), p (t = 0)] (510) where, the voltagedependence of the functional is eliminated. Based on the functional I[Pe, ,Pn j, all steady state properties can be expressed as functionals of Pe, and p,. 5.2.3 The Stationary Principle and the Single Particle Picture For the timedependent system, we can define the action A as, A = ~(t)i IH(t) V(t))dt (511) where y,(t) is the many body wave function. If y(t) is normalizable, then given the initial state, y(oo), the stationary requirement of the action, 6A =0 (512) 9(0 is completely equivalent to equation (53), where 38(t) is an arbitrary change of the wave function under the conditions that 38/(0) = 0 and the normalization is unchanged. There is no difficulty to normalize the wave function in the closed system, but it seems to be impossible to normalize the wave function for the open system which has the geometry shown in Fig (11). Thus, the Equation (512) can not be directly applied. However, for the case of timedependent DFT, which means that the initial state is a stationary ground state, the action is the functional of the electron density A = A[p]. As long as the wave function y(t) yields the correct electron density, we still can use the Equation (512) to calculate the action and the stationary requirement becomes, 9A = 0 (513) 8p(t) with the conditions that 8p(0) = 0 and 8p doesn't change the neutrality condition of two leads and the scattering region. This procedure can be generalized to the nonequilibrium initial state case. According to the previous subsection, for the open system under a finite bias, the action can be written as A[peq, np, and the stationary principle becomes, 6A 6A =0 (514) pe, (t) 8pn (0 where the arbitrary change of electron densities should not change the initial conditions nor the neutrality condition of the system. For the adiabatic timedependent process, the action reduces to 0 A= E(t)dt. (515) Then, the Equation (514) becomes, = 5E0 (516) 3p,q (t) gp, (t) where, E is the total energy of the system at time t. Both the ground state and the timedependent DFT work in the framework of a singleparticle picture. We follow the same procedure. Consider an independent particle system, which has a same geometry as shown in Fig. 11 and whose single particle orbitals can be written as #,f, where m denotes all good quantum numbers other than the energy At the initial time, since the two leads are wellseparated with respect to the scattering region, fm, consists of three kinds of states, 0,,, and ;b, which are localized in the left lead, the right lead, and the scattering region respectively. Those three kinds of states obey the equilibrium distribution functions of three regions respectively. When the temperature is zero, those distribution functions are step functions. The adiabatic timedependent procedure described in previous sections will not change the occupation number of each state when it evolves with the time. Further assume that at the initial time, 0, ,,0, and b ,, are perfectly normalized such that the neutrality conditions are satisfied in all three regions. Since the system is noninteracting, both Green's functions G (E) and G' () are diagonalized by m,~ It is straightforward that equation (57) becomes, Peq= Z ,e 2 m2 (517) P, = fm( )W,2R m,/L >E>/PR where f, is the occupation number. Since the occupation number of each state does not change with time and/ u > R,, p,, is, P, = 2 (518) Assume that Peq and p,, yielded by equaton (517) and (518) are the same as those of the true interacting system under study. Then, similar to the equilibrium state DFT, we write the total energy asE = T + E where T is the singleparticle kinetic energy, T = T +T Teq Z (O V2 m). (519) m,E<,uR T,= Z (, Iv20) l,/L >E>/R 2 The potential energy E, can be written as, E= rV, (r)p(E) +Ec, +E E 1 J 1 (520) 2: j 2 where EcoU, is the classical electronelectron repulsion energy and Exc is the exchange correlation energy. E, is the functional of peq and p, We can define two effective potentials based on E, as, Jeq 8EV Vef  Peq .(521) Ve" = E ef Pn The Equation (516) can be written in terms of those orbitals as, (E V2 ), R = 0 2 (522) 1 (E + V2 2 V,,) m,/i.>>>R = 0 There are two boundary conditions for the system under study. The open boundary condition incorporated with equation (522) yields the scattering states and the bound state boundary condition yields bound states. Unlike equilibrium DFT, we have two different effective potentials instead of one. The nonequilibrium electrons see a different mean field potential than the equilibrium electrons do. The selfconsistent procedure makes two electron densities Peq and p, to converge respectively, and the current is carried by the nonequilibrium electrons. 5.3 Exchange energy for uniform electron gas The exchangecorrelation energy can be written as the sum of the exchange part and the correlation part, Ex = Ex + Ec. Both Ex and Ec are functionals of Peq and p,. They are very complicated for a real, inhomogeneous system. In this section, as a case study, we calculate the exchange energy for the uniform electron gas under a finite bias, which represents a perfect conductor connected to two reservoirs infinitely far away. This study serves as a test for the theory described in the previous section. For the uniform system, if we take periodic boundary condition as q(x + 1) = O(x), then the singleparticle orbitals can be written as, (k, k, k) = e (kx+kk ) (523) 2V2 /3/2 W12 where kxy = nx,y, with nx,y, are integers. Set the z axis to go from left to right, which is the direction of the current. There are two kinds of states for this system: the rightgoing states (k > 0) and the leftgoing states (k < 0). From equations (517) and (518), we see that the occupation of those states is: belowke, all states are occupied and between k, and k,, only rightgoing states are occupied. Define two quantities k, and k, h2k2 h2k2 by = /R, = L. Two electron densities P,, (which comes from all states below 2m 2m ke) and p, (which comes from all rightgoing states above ke), can be related to k, and k3 1 k,. After some algebra, we have peq = and p,,6 (k k). We can see more 37T 6;T clearly from this case that the total electron density p = Pe + p, alone is not enough to describe the system: two systems with the same p have different electronic structures if they have different set of Peg and p,. The exchange energy Ex can be evaluated by the equation [30] E = 1 ri(\,r) cI2clrl (524) 4 r ll where (ri,2 ,) is the density matrix and its diagonal part is the electron density. It can be calculated as p(r, r) = 2 Yer'12 where the summation is over all occupied k For a large number of occupied states, the summation can be replaced by an integral, S ke 2T I kt ,T/2 2, p(, i) = k2dk SinOdO fdke'12 + k 2dk f SinOdO [de 12 (525) 0 0 0 ke 0 0 where the first integral is the contribution of all states below ke, which is an equilibrium state integral and can be evaluated analytically. We denote this term as Pe. The second integral is the contribution of occupied states between k, and k,, which can be denoted as /,. For the nonequilibrium integral, since only states with positive kZ are occupied, the variable 0 of the second integral goes from 0 tor /2 We introduce two variables R = (tr + 2)/2 and S = '2. Then the density matrix is only a function of s and the exchange energy can be written as, E = dRK 1 1 (526) K= p(s) 2 4 s where K is a constant since the system is uniform. For a given electron density p, when the bias voltage increases, p,, increases and Peq decreases until at a certain bias voltage Vmax, Pe becomes zero and p,, equals to p. Vmax defines the biggest difference between /, and /R we can have for a uniform system with a given p at zero temperature. Given the electron density p and the bias voltage L /R we can numerically evaluate p andK. Fig. 51 shows the ratio K/Ko as a function of total electron density and bias voltage. Ko denotes the term calculated at zero bias. From the figure, we clearly see the voltagedependence of the exchange energy for a given electron density p. This supports the analysis of the previous section: the electron density alone is not enough to describe the system. Also, the p dependence of the exchange energy at a given bias voltage, which will determine the form of the density functional, varies with the bias voltage. It means that we will have difficulty in finding a universal functional for all bias voltages, which agrees with the previous section too. 1.012 1.008 / / o \ \ p= 3377 1.004 \ p=17.29 1.000 p=7.30 0.996 =2.16 0.992 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Bias Voltage (Volts) Figure 51. Exchange energy at different bias voltages for different p (Electrons/m3) We notice that although the voltagedependence of the exchange energy varies with the electron density, the behavior is similar for all electron densities: as the voltage increases, the exchange energy increases until it reaches a maximum and then decreases. In order to understand this behavior better, we rewrite the Eq. (526) as, K =Ke +K +K Keq f Ieq 1 :1 2 (527) K = 4 s 2 4 s The first term Keq is the contribution of the exchange of equilibrium electrons. The second term Kn comes from the nonequilibrium electrons and the third one denotes the contribution of exchanging between equilibrium and nonequilibrium electrons. The voltagedependence of these three terms for p = 7.29 Electrons/m3 is shown in Fig. 52. At bias zero, since p, is zero, there are no nonequilibrium electrons. Hence K, and Keqn are zero. As the bias voltage increases, peq decreases, which leads to the decrease ofK, and p, increases, which leads to the increase of both K, and Keqn The term Keq n reaches its maximum at a certain bias voltage and then decreases. At the maximum bias voltage, p, = p and there are no equilibrium electrons, then both Keq and Keqn are zero at this bias voltage. According to the previous section, the exchange energy can be expressed as a functional of p and p,. We write K as, K = Ko [p1]+ f[ p, p (528) where, Ko [o] is the one for the equilibrium system with the electron density p; this term is well known. We plot the numerical results of K / K as the function of p, / p for different electron densities in Fig. 53. 80 12 i p12 =7.30 0 10 E0 o x 8 E" 0 0 ,, \ / 2 o6 OO Lo 2 O m m O) 3.0 =7.30 m m m m m 1.5\  1.0 0 (b) 3.0 p=7.30 2.5  2.00 electrons and the exchange between nonequilibrium and equilibrium Sa a 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Bias Voltage (Volts) (a) 3 .5  2.0Figure 52. Split of the exchange energy: (a) the exchange between nonequilibrium 0.5electrons and (b) the exchange between equilibrium electrons. electrons and (b) the exchange between equilibrium electrons. The numerical results show that we have a universal functional f[ for all I P electron densities. These data can be fit to the function, f PJ +b (529) where a = 0.0497 and b = 0.0436. 1.012 1.008 1.004 1.000 0.996 0.992 0.2 0.4 0.6 0.8 p' p Figure 53. Exchange energy as the functional of p, / p for different electron densities. Then the exchange energy can be written as, E=C 4/3 + P (530) where, C, = 0.7386. Equation (530) gives the exchange energy functional of the uniform system for all bias voltages. It clearly shows the dependence on the nonequilibrium electron density as well as the dependence on the total electron density, which support the analysis of the mp=7.30 0 p=2.16 + p=17.29 x p=33.77 previous section. Whether or not this functional can be generalized to the inhomogeneous system has to be resolved in a future study. 5.4 Conclusion In conclusion, we have shown that the transport properties of an open system under a finite bias can not be determined by the total electron density alone. We introduce two electron densities pe and p, to describe such a system and show that the non equilibrium electrons experience different effective potential from the equilibrium electrons. As a test, we study the exchange energy of the uniform electron gas under a finite bias. The numerical results clearly show the dependence to the nonequilibrium electron density, which support the theory. We fit the numerical data to the exchange energy functional. The applicability of the functional to the inhomogeneous systems needs to be addressed in future studies. CHAPTER 6 SUMMARY AND CONCLUSIONS We introduce the main ideas of first principles method for open system under a finite bias. The method has been successfully applied to study the transport properties at nano scale, for example, molecular electronics. The achievements of the method can be summarized as follows. The shape of calculated IV curve can be compared to the experiments and peaks in the curve can be explained by tunneling through molecular orbitals. We apply the stateofart approach which combines the nonequilibrium Green's functions technique and DFT to Au/S/Azobenzene/S/Au molecular junction. The results show that the azobenzene molecule is a good candidate of lightdriven molecular switch, which is important for application of nanoelectronics. The IV curve of cis molecular shows a negative differential conductance at a certain bias, which is due to the suppression of the tunneling through HOMO due to bias voltage. We also generate the method to include spin and implement the method in a wave function representation. The method is used to study the TMR through a layered Fe/FeO/MgO/Fe magnetic junction. The results suggest a biasenhanced TMR which is rarely seen in previous experiments. The enhancement of TMR due to bias voltage is explained by the effects of tunneling through the surface resonance states. Our analysis suggests that this phenomenon should be able to seen in other systems. In the final part of the thesis, we discussed the problem in the current approach and suggest that in order to correctly describe the nonequilibrium effects introduced by the 84 external bias, an additional parameter other than the total electron density should be introduced into the equilibrium DFT. For a case study, we calculate the exchange energy for an uniform electron gas and clearly see the voltagedependence of exchange energy. By introducing 'nonequilibrium electron density', we show that the voltagedependence of the functional can eliminated. This work is supported by the Alumni Fellowship of University of Florida and the US Department of Energy under contract DEFG0202ER45995. APPENDIX A EXPANSION OF LESSER GREEN'S FUNCTION The nonequilibrium Green's functions are defined as G(1,1') = i( )(1 (Al) where, the shorthand notation 1 (x, t,) is used. The timeorder operator Tc is defined on two time branches, + branch, which runs from t = co to t = +oc, and branch, which runs in the opposite direction. With two time labels, tj and t,, which can be located on either of two time branches, we have four kinds of nonequilibrium Green's functions: G(1,) t, t, e+ G'(1,1') tl e,t, e+ G(J1,') = (A2) G'(1,1') tj +, t,e G' (1,1') ti, ti,  These four Green's functions can be denoted as the timeordered Green's function Gt, the greater Green's function G', the lesser Green's function G', and the antitime ordered Green's function G'. For the Hamiltonian which is given by equation (21), a particular lesser Green's function Gk defined as G (, (A3) 86 where, the indices k, a denote channels of the left lead and n denotes the basis states of the device region. The angle brackets indicate the ensemble average. The time of the operator d+ lies on the branch and the time of the operator Cka lies on the + branch. tka tV + tk. v 9kta,ka ka,m mGn tn 1 tn tkc tv 'LW g ka kc kcxm mn m Figure A1. Expansion of lesser Green's function Gkx,,. The time of the perturbation can be either on + branch or on branch. The coupling term Hoo, is taken to be the perturbation V Then, the Dyson equation G = G, + GVG can be used to expand the Green's functions, where the term Go is the Green's function without the perturbation. Inserting the basis states into the Dyson equation, we have, Gk,,, = (ka G n)= (ka G n) + (ka GVG n) = gkz, mG,,,, (A4) where, the term (ka GG, n) is zero because electrons can not propagate from the left lead to the device region without the coupling term and the term gk,,k is the bare Green's function of the left lead. The equation (A4) can be used for both equilibrium and non equilibrium Green's functions. For nonequilibrium case, since the time of the perturbation can be either on + branch or branch, the equation (A4) gives two terms. 87 As an example, the expansion of the lesser Green's function Gk contains two terms as shown in Fig. A1. Thus, Gk can be written as Gkn = 7Vkam kakaGn gakaGmn. (A5) m APPENDIX B DETAILS OF NONEQUILIBRIUM SELFCONSISTENT PROCEDURE As mentioned in Chapter 3, the device region includes one unit cell of the left lead (denoted as L), the molecule (denoted as M), and one unit cell of the right lead (denoted as R). The first step of the calculation is to compute selfenergies of two leads according to equation (32). Since the tightbinding model is used to describe two leads, only the part L (R) of the device region is affected by the left (right) semiinfinite lead. Thus the dimension of the selfenergy of the left (right) lead YL (R ) is the number of basis states of one unit cell of left (right) lead. In our calculation, the Gaussian basis set is used. Due to screening, selfenergies of two leads are not affected by the device region, so they don't change in the selfconsistent procedure. The selfconsistent procedure starts from an initial electron density p, (r) in the device region. Given the electron density, the Hamiltonian can be calculated as H 1v (') H 2 2 +'+ V (r)+ + Vbl]+ (r), (Bl) 2 Pr where, the last term Vb,, (r) is the linear voltage drop, Vb,K, (r) = R y The current d is flowing in y direction and d measures the distance between two leads. The retarded and advanced Green's function can be achieved by Gr', = (B2) SS t t where, S is the overlap matrix of Gaussian basis states and H is the Hamiltonian matrix. In order to see it more clearly, equation (B2) can be expanded as K+SLL HLL HL +SML HML E+SRL HRL G' = E '+S HLM ES Ha ~+ Sm H H (B3) E + LR HLR +S H E+S H~ RR R where, indices L, M, and R denote three parts of the device region respectively. Then, the lesser Green's function can be computed according to equation (211) and the density matrix is determined by jP, = IG; (8)dE, (B4) where, indices i, j denote Gaussian basis states. The electron density can be directly related to the density matrix as p()p) =( = )= (r i)(i p j)(j 6) = 01, ()p,,b0(), (B5) 1,j 1,j where, 0, (i) is the Gaussian basis in real space. The electric static potential is defined as S() () + Vb, () (B6) An interesting quantity is the difference between the electric static potential under a finite bias and the electric static potential of zero bias, Ve , AV,, = V ,, V, = Vb, (i) + tep(c') (B7) which shows the effects of the voltage bias on the electric static potential. which shows the effects of the voltage bias on the electric static potential. 