Citation
Transport properties at nano scales via first principles studies

Material Information

Title:
Transport properties at nano scales via first principles studies
Creator:
Zhang, Chun ( Dissertant )
Cheng, Hai Ping ( Thesis advisor )
Hebard, Arthur F. ( Reviewer )
Herschfield, Selman ( Reviewer )
Roitberg, Adrian ( Reviewer )
Trickney, Samuel B. ( Reviewer )
Place of Publication:
Gainesville, Fla.
Publisher:
University of Florida
Publication Date:
Copyright Date:
2004
Language:
English

Subjects

Subjects / Keywords:
Atoms ( jstor )
Electric current ( jstor )
Electric potential ( jstor )
Electrodes ( jstor )
Electron density ( jstor )
Electrons ( jstor )
Greens function ( jstor )
Lead ( jstor )
Molecules ( jstor )
Statistical bias ( jstor )
Physics thesis, Ph.D
Dissertations, Academic -- UF -- Physics

Notes

Abstract:
There are two main difficulties for the first principles study of transport properties at the nano scale. The first is that many-body 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 non-equilibrium 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 Non-equilibrium 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' many-body effects and NGFT naturally includes the non-equilibrium effects. The new approach uses NGFT to treat the non-periodic 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 light-driven molecular switch. The switch, made of a single molecule, is one of the most important elements of nano-electronics. However, most proposed molecular switches are driven either by an external bias voltage or by STM manipulation, neither of which is ideal for nano-scale circuits. The switch we designed has a high on-off conductance ratio and more importantly, can be driven by photons. In following chapters, we generalize the method to the spin-dependent case and apply it to a magnetic layered structure. We implemented the method within the framework of the Layer Korringa-Kohn-Rostoker (LKKR) approach, which is particularly well-adapted to the layered structure and found a bias-enhanced tunneling magneto-resistance (TMR) for the Fe/FeO/MgO/Fe junction. Our results are important not only for application, but also for understanding of the voltage-dependence 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 impurity-free layered structure, a different behavior of TMR may occur.
Subject:
density, light, magnetoresistance, molecular, nanoelectronics, nonequilibrium, quantum, spin, transport
General Note:
Title from title page of source document.
General Note:
Document formatted into pages; contains 105 pages.
General Note:
Includes vita.
Thesis:
Thesis (Ph.D.)--University of Florida, 2004.
Bibliography:
Includes bibliographical references.
General Note:
Text (Electronic thesis) in PDF format.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright Zhang, Chun. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:
12/18/2004

Downloads

This item has the following downloads:


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. Hai-Ping 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 long-term 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 Landauer-type Formula ..................................................4
2.2 The Combination of DFT and NGFT for An Open System ...................................8
2.2.1 Self-Consistent Procedure for Ground State of Closed System ...................8
2.2.2 Self-consistent Procedure for Open System ........ ........... .................. 10

3 COHERENT ELECTRON TRANSPORT THROUGH AN AZOBENZENE
MOLECULE: A LIGHT-DRIVEN MOLECULAR SWITCH............... ...............12

3.1 Introduction ................................................. ........ .. ....................... 12
3.2 Tight-binding M odel and Self-energies...................... ..... ................. 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 Non-equilibrium Green's Function Formalism in Real Space ...........................34
4.3 Self-consistency 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 NON-EQUILIBRIUM APPROACH FOR OPEN
SYSTEMS BASED ON DFT .............. .............. .... ............... 63

5.1 Introduction .................. ...... ............ ....................... ................ .....63
5.2 Hohenberg-Kohn 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 NON-EQUILIBRIUM SELF-CONSISTENT 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

1-1. Partition of the open system connected to two metal electrodes ............................2

3-1. Physical models of an azobenzene switch in contact with two gold leads: trans
isomer (upper) and cis isomer (lower). ........................................ ............... 14

3-2. Tight-binding model of an open system connected to two semi-infinite leads. Indices
i, j denote nodes of device region and a, b are the nearest node of left and right
lead. Only the nearest-neighbor interaction is considered. .....................................15

3-3. A perfect lead that consists of two semi-infinite leads and a hopping term. Indices m
and i belong to the left semi-infinite lead and i+1, j belong to the right semi-infinite
lead ................................................................................. 16

3-4. 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

3-5. 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

3-6. Local DOS at the Fermi energy projected on seven sites in the switch. From 1-7, or
from left to right, Au lead, S atom, benzene ring, double-bonded N2, another
benzene ring, S atom and Au lead ............................... .....................22

3-7. Non-equilibrium self-consistent procedure. The non-equilibrium Green's function is
calculated according to equation (2-11). ....................................... ............... 24

3-8. The difference of electronstatic potential between bias one volt and bias zero for cis
configuration. The x-z plane is parallel to the N2 bond, cutting through the
m idpoint of A u unit cell. ................................................................. .....................25

3-9. The transmission probability as a function of single-particle energy and the external
bias voltage for cis configuration. ........................................ ......................... 26

3-10. 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









3-11. The differential conductance as a function of bias voltage for both trans and cis
configurations............................................................................................. 28

3-12. The transmission as a function of energy at 0.64, 0.72 and 0.8 Volts for cis
configuration ...................................................... ................. 2 8

4-1. Shifts in the density of states in each atomic layer at a bias voltage of 0.54 V..........44

4-2. 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

4-3. 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

4-4. 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

4-5. 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

4-6. 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

4-7. 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

4-8. 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

4-9. 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

4-10. 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









4-11. 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

4-12. 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

4-13. 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

5-1. Exchange energy at different bias voltages for different p (Electrons/m3) .............78

5-2. Split of the exchange energy: (a) the exchange between non-equilibrium electrons
and the exchange between non-equilibrium and equilibrium electrons and (b) the
exchange between equilibrium electrons. ..................................... ............... 80

5-3. Exchange energy as the functional of p/ p for different electron densities ..........81

A-1. 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: Hai-Ping 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 many-body 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 non-equilibrium 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 Non-equilibrium 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' many-body effects and NGFT naturally includes the non-

equilibrium effects. The new approach uses NGFT to treat the non-periodic 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

light-driven molecular switch. The switch, made of a single molecule, is one of the most









important elements of nano-electronics. However, most proposed molecular switches are

driven either by an external bias voltage or by STM manipulation, neither of which is

ideal for nano-scale circuits. The switch we designed has a high on-off conductance ratio

and more importantly, can be driven by photons. In following chapters, we generalize the

method to the spin-dependent case and apply it to a magnetic layered structure. We

implemented the method within the framework of the Layer Korringa-Kohn-Rostoker

(LKKR) approach, which is particularly well-adapted to the layered structure and found a

bias-enhanced tunneling magneto-resistance (TMR) for the Fe/FeO/MgO/Fe junction.

Our results are important not only for application, but also for understanding of the

voltage-dependence 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 impurity-free 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. [1-5]. 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 time-dependent

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 steady-state 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 steady-state 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 effective-mass assumption. The Boltzman equation is based on a

probabilistic description for the evolution of a single-particle distribution f(ri, p, t)

which is a function of classical variables r and p However, when the device region is at

nano-scale, 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 non-equilibrium

Green's functions technique (NGFT) and Density Functional Theory (DFT) has proven to

be successful [6-22]. The method can be summarized in two steps. First, a non-

equilibrium self-consistent procedure based on DFT and NGFT is used to determine the

single-particle effective potential. Second, a Landauer-type formula is used to calculate

the current, which relates the current to the integral of the single-particle transmission

probability over energies.


Current I




Left lead Device Region Right lead



Figure 1-1. Partition of the open system connected to two metal electrodes.

The approach combining NGFT and DFT is based on the single-particle 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 two-step 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 light-driven

molecular switch.

In Chapter 4, we generalize the method to spin-dependent tunneling junctions with

application to a layered magnetic Fe/FeO/MgO/Fe junction. The calculation suggests a

bias-enhanced tunneling magneto-resistance (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 two-step approach.

For the first step, in which the effective single-particle potential is determined, DFT is the

biggest problem. The Hohenberg-Kohn 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

non-equilibrium effects into DFT. For the second step, in which a Landauer-type formula

is used to calculate the current, the problem is that the Kohn-Sham 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 Landauer-type 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 non-interacting, 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 (2-1)

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 non-interacting, we have

HLR = kak,a
k,aeL,R

en = nt({d, {d) (2-2)

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 single-particle 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 (2-3)


where NL = 'c-ck, 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
(2-4)
e knka Vka,nG
k,aGL
n

where the current is expressed in terms of two non-equilibrium '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 (2-5)
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 non-interacting 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 (2-5) into equation (2-4) and doing a Fourier transform gives


J, = 2e d Re{ VknVk m k ()Gm (0) + gK~ ,ka ()Gn( (o)]} (2-6)
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 (2-7)
gka-,ka ka i


where f(ska) is the Fermi-Dirac distribution function. Insertion of the equation (2-7)

into equation (2-6) gives the current as


J = ie Tr{L [fL (Gr G) + G]} (2-8)
f 2;r

where the line-width function FL is defined as

[F (k )]mn = 2r ( pa (c )V. (,k )Vl (k) (2-9)
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]} (2-10)
2 f 2;r

where the current is related to the local properties of the device region.

In equation (2-10), 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 non-interacting. For the free-electron model, we have [15, 20,

25]

G' = Gr''G
Y< = i[F/f, +FRR (2-11)

and then the current can be written as

J= elde[f fR ]Tr{GTRGTr } (2-12)

which is a Landauer-type formula.

Since the interaction in both leads and the device region usually plays an important

role in the tunneling process at nano scales, equation (2-12), 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 single-particle picture. Within

present-day approximations, it is able to take into account 'most' many-body effects.

Moreover, due to the success of DFT in the electronic structure calculation of molecules,

it is natural to combine DFT and the Landauer-type formula to study molecular

electronics.









2.2 The Combination of DFT and NGFT for An Open System

2.2.1 Self-Consistent 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 [26-30]. The main ideas of DFT are outlined in the

following paragraphs.

The basis of DFT is the Hohenberg-Kohn 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 (2-13)

where the term Vex, represents the entire external potential (including the ion-electron

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 ground-state electron density minimizes the energy functional E[p]. The only thing

used in the proof of Hohenberg-Kohn Theorem is the minimum energy principle of

ground state.

On the basis of the Hohenberg-Kohn theorem, a self-consistent procedure based on

the single-particle picture can be constructed. Assume that there is a non-interacting

system which yields the same electron density as the interacting system of interest, then

the electron density can be written in terms of the single-particle functions

p(/) = ~ (F) 2 (2-14)
1









which are usually called Kohn-Sham orbitals. The kinetic energy of the non-interacting

system can be calculated in terms of the Kohn-Sham orbitals (Hartree atomic units)

V2
T/[pT]= -- ). (2-15)

Then the exchange-correlation energy can be defined as

Exc [p] = F[p] T,[p] EH [p] (2-16)

where the term EH [p] is the classical Coulomb energy. The exact form of Ex[p] is

unknown. However many kinds of exchange-correlation functionals, most based on the

local density approximation, have been proposed and successfully used. By minimizing

the energy functional, we get the Kohn-Sham equations

V2 + VH [p](r) + Ve, () + V [P](ir), (i)= g () (2-15)


where the first term VH is the well-known classical Coulomb potential and the term Vxe


is the exchange-correlation potential defined by [-- The effective Kohn-Sham
9p

potential is defined as

VeY [p]() = VH [p](r) + Ve, (r) + VX [p]() (2-16)

Now, a self-consistent procedure (Kohn-Sham procedure) can be completed as follows.

Given the electron density p, the effective Kohn-Sham potential Vff [p] can be calculated

according to equation (2-16). Knowing VJ. [p], we can get the Kohn-Sham orbitals V, by

solving equation (2-15) and then the new electron density can be calculated according to

equation (2-14). Repeat the procedure until the self-consistency is achieved. The

Hohenberg-Kohn theorem guarantees that the Kohn-Sham 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 self-consistent

procedure based on the single-particle picture.

2.2.2 Self-consistent Procedure for Open System

For an open system with the geometry shown in Fig. 1-1, 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 Kohn-Sham 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 Kohn-Sham 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 self-consistently 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 (2-17)
2n i









For a single-particle Hamiltonian, the 'lesser' Green's function can be calculated by the

equation (2-11), 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 (2-18)
E C y r(a) r

In equation (2-18), the term HC,, denotes the effective single-particle Hamiltonian and

the two terms EL ,R are self-energies 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 self-consistent 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 (2-17).

After self-consistency is achieved, the Landauer-type formula (2-12) can be used to

calculate the current in the system.

The non-equilibrium self-consistent procedure mentioned above is implemented in

a Green's functions representation which employs a multi-band tight-binding 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 LIGHT-DRIVEN MOLECULAR SWITCH

3.1 Introduction

Electronic devices soon will be constructed of molecules or have features of

molecular size. The results of recent experimental [1-5] and theoretical studies [6-20]

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 [31-35]. 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, light-driven

processes, are therefore highly desirable. The photosensitive azobenzene molecule

discussed here is one such candidate for an opto-electronic device.

The azobenzene molecule has two stable conformations [36-38], 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 opto-mechanical 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, Au-S-

azobenzene-S-Au. 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 semi-infinite gold leads. Note that the name of the modified

molecule is technically bis-(4-methanethiol-pheyl) 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. 3-1 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 non-equilibrium case. In both cases, the localized

basis set and the tight-binding model is employed.








k=365 nm X=420 nm







Figure 3-1. Physical models of an azobenzene switch in contact with two gold leads:
trans isomer (upper) and cis isomer (lower).

3.2 Tight-binding Model and Self-energies

As mentioned in Chapter 2, the key issue is to evaluate the self energies of two

semi-inifinite 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 self-energy term represents the effects of the semi-infinite 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 tight-binding model is a natural choice. Fig. 3-2

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 3-2. Tight-binding model of an open system connected to two semi-infinite leads.
Indices i, j denote nodes of device region and a, b are the nearest node of left
and right lead. Only the nearest-neighbor 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, (3-1)

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 (3-2)
R = h gbbh

where, EL, 2R are self-energies of the left and right lead respectively. Two hopping terms

are projections of the effective Kohn-Sham Hamiltonian to a localized basis of each node

(hopping terms are matrices also). The only problem remaining is two so-called surface

Green's functions ga and gbb

Since a perfect lead with periodic boundary conditions consists of two semi-infinite

leads and a hoping term within the framework of tight-binding as shown in Fig. 3-3, it is

possible to construct the Green's functions of semi-infinite 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 3-3. A perfect lead that consists of two semi-infinite leads and a hopping term.
Indices m and i belong to the left semi-infinite lead and i+1, j belong to the
right semi-infinite 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 (3-3)

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,+, (3-4)

Similarly,

G,,, = G y G, = g+ ,zh G, ,. (3-5)

In order to solve g from equation (3-4) and (3-5), 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 tight-binding Hamiltonian is

H = hoa,, + C haa,, + h a+, am, (3-6)
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, (3-7)

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 (3-7) can be written as

Az -h h- z 2=, =0, (3-8)

where, M+ = E+ ho and h, = h E It is to be noted that not all eigenvalues A, of

equation (3-8) for a given energy represent the Bloch state. Because equation (3-8) 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 (3-9)
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 =
(3-10)
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). (3-11)









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. (3-12)

Comparing with equation (3-10), we see that, for i < j, A = a Introducing

G+ = AG+, into equation (3-5), we have

g, = a (h) (3-13)

Similarly, fori > j, we have

g ,+, = p+ (h+ ) (3-14)

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 linear-response 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 ) (3-9)
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, STO-3G

Gaussian basis [28] for all other atoms (C, S, N, H). The pbe0 exchange-correlation

functionals [26] are used. For the wire system depicted in Fig. 3-1, the system is, in

principle, infinite, and periodic boundary conditions are not appropriate. However, due to










screening, the Kohn-Sham potentials in the two semi-infinite 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 3-4. 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 self-energy 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] (3-10)



which is the kernel of the integral in equation (2-12). 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 3-5. 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. 3-4. Unlike the

DOS of an isolated molecule, Fig. 3-4 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. 3-5, 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. 3-5 indicates clearly that the broad

background in the DOS derives primarily from the Au layers.

Figure 3-5 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 zero-bias 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. 3-4.

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.3-4, 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 zero-bias conductance for the trans


structure is 0.21 x 104 4-1 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 3-6. Local DOS at the Fermi energy projected on seven sites in the switch. From
1-7, or from left to right, Au lead, S atom, benzene ring, double-bonded N2,
another benzene ring, S atom and Au lead.

To understand the tunneling process at the Fermi energy, Fig. 3-6 depicts the local

DOS at the Fermi energy, Fig. 3-6 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

double-bonded 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 non-equilibrium self-consistent 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 zero-bias case. We neglect the effects of the so-called current-induced

force by making this assumption.

1 G
2; i





Vff = V[p]

Figure 3-7. Non-equilibrium self-consistent procedure. The non-equilibrium Green's
function is calculated according to equation (2-11).

The self-consistent procedure is depicted in Fig. 3-7. 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 non-equilibrium Green's

function G' can be calculated according to equation (2-11), where the retarded and

advanced Green's functions are calculated in the same way as in the previous section.

The details of the self-consistent procedure can be found in Appendix B.

One of the major accomplishments of the non-equilibrium first principles

calculation is that the electrostatic potential due to the external bias can be computed self-

consistently. Fig. 3-8 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 x-z 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







0-1



S(Angstros 10


Figure 3-8. The difference of electronstatic potential between bias one volt and bias zero
for cis configuration. The x-z plane is parallel to the N2 bond, cutting through
the midpoint of Au unit cell.

The current through the system is calculated by equation (2-12), where the current

is expressed as the integral of the single-particle 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 3-9, 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)).


(3-13)









Transmission as a function of energy and voltage for cis

1.0


HOMO
c / LUMO







0.0Ie 7 0\ \


nergyV) -4 00 ,9


Figure 3-9. The transmission probability as a function of single-particle 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 Kohn-Sham

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.

I-V curves through the trans and cis configuration are shown in Fig. 3-10. The

differential conductance, defined as the derivative of the current with respect to the bias

voltage, is shown in Fig. 3-11. 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 3-10. 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. 3-10. 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. 3-11. This can be explained by consideration of

Fig. 3-12, 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 (3-1), 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. 3-10 clearly shows such an example, where as the bias voltage
































0.0 0.2 0.4 0.6 0.8 1.0
Volts


Figure 3-11. 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 3-12. 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. 3-11.














CHAPTER 4
QUANTUM TUNNELING THROUGH MAGNETIC JUNCTIONS

4.1 Introduction

Quantum transport in open systems far from equilibrium is an important and

challenging problem [1-5] 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 [6-20] with various

degrees of success. The preceding chapter is an example. The problem of spin-dependent

tunneling in magnetic heterostructures, however, presents a significantly greater

challenge than the spin-independent case to the theory of non-equilibrium transport

because one must simultaneously deal with the complexities of non-equilibrium transport

and complex electronic structures along with magnetism. Spin-dependent tunneling

junctions are considered to be the most likely candidate for next-generation devices for

spintronics. Compared to the trilayer or multilayer giant magnetoresistance (GMR)

systems, magnetic tunneling junctions tend to have a relatively higher magnetoresistance

ratio [45-48], which is an important factor in device designs. The magnetoresistance ratio

is defined as (Gp-Ga)/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 anti-parallel. One of the









significant differences between a GMR system and a spin-dependent 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 first-principles studies [50-54] 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 linear-response 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 [51-52]. 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 voltage-dependence of TMR in

terms of the impurity scattering at the interface. They suggested that with a high quality

barrier, a different voltage-dependence of TMR may occur. Some groups [55-57] 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, first-principles calculation [6-20] of electronic structure and tunneling

conductance under a finite bias voltage has been successfully applied to molecular

electronics. As mentioned in Chapter 1, the first-principles approach can be built on

either non-equilibrium 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 Korringa-Kohn-Rostoker (KKR) approach [61], which was used

successfully in calculating linear-response transport properties of spin-dependent

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 two-terminal conductance within the linear-response 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 [62-65]. The

inclusion of a single atomic layer of FeO on one of the interfaces is important. An earlier









linear-response 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 first-principles 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 exchange-correlation 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 exchange-correlation potential.

In the following section, we will show that our approach is equivalent to the

Keldysh Green's function approach within the single-particle 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 Non-equilibrium 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)) (4-1)

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 (4-2)
G'(,';E) -G'(F, ';t)e-EtdAt

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. (4-3)
Z.7U> '









The problem remains to evaluate G' (r,'; E). For non-interacting 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. (4-4)


In the framework of mean-field theory, these two operators can be expanded as single-

electron operators

+ (ir, t)= Y0(r, t)C +
n (4-5)
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 (4-5) into equation (4-4), we have

G< (, ',E)= i 9(E- E) J Cn )CC n (r')o (r)
(4-6)
= 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 (4-6) 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 1-1. 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 Fermi-Dirac 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 (4-6) 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). (4-7)
1 m

We compare this result with that obtained in Chapter 2. Rewrite equation (2-11) as

G' = fG'FLG" +fGrFRGa. (4-8)

Both equations are valid within the single particle picture. However, equation (4-7) is

more general than the equation (4-8), 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 (4-7) applies to the entire space whereas

equation (4-8) can be applied only within the device region.









Substituting equation (4-7) into equation (4-3), we obtain the scattering part of the

electron density

p, (r) = f (E k) 2 + f (E R) 2. (4-9)
/ 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 single-particle functions of bound states. Accordingly, the electron

density also has an extra contribution from


p Zf (E)b 2 (4-10)
b

The determination of fb is not unique and is still being investigated by researchers.

4.3 Self-consistency and Tunneling Current for A TMR Junction

In this section we outline our approach for non-equilibrium self-consistent

calculations of electronic structure of a TMR junction.

Although the approach applies generally to any lead-device-lead ensemble, the

particular Fe/FeO/MgO/Fe tunneling junction in our calculation has a layered structure

with a two-dimensional periodicity along each layer. This approach allows a two-

dimensional Fourier transform that yields a quasi one-dimensional problem for each

reciprocal vector k//. The left and right leads are semi-infinite 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 right-going states

(+ and left-going 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 (4-11)

to (0' Right lead
m

and

( +-t ,t' ,,,(+', Right lead

S= ( device region (4-12)
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 (4-9) and (4-10) we obtain


p(r) = dEp (r) + dEp (r) +pb(r) (4-13)


where


p =(r) YZ (E- E,) ,
2 (4-14)
p (r) = Z (E Em) (
m

These equations hold for all regions including the device region and both leads. We can

rewrite equation (4-13) as


p(r) = p (r) + dEp (r) (4-15)
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) (4-16)


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 (4-15) 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) (4-17)
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) (4-18)
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 self-consistent procedure for the non-equilibrium 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 semi-infinite (lead) parts are

calculated using the bulk self-consistent 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 spin-dependent calculation is done, we write the electron

density in equation (4-14) in terms of spin explicitly


PE,(r)= Z8(E- Ek)
2 (4-19)
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) (4-20)
PlR

where the equilibrium part is


pq, (r) = JdE[p (r) + PE (r)]+ p, (r). (4-21)


The last term means only those bound states which have spin c are counted.

The equation (4-21) 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

spin-density functional to calculate the effective potential.

Similarly, the current density for each spin is given by a Landauer-type formula


J = eJdE e ti m (4-22)
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 (4-20) can be evaluated using the standard

contour integration technique for Green's function methods. The non-equilibrium term on

the other hand, is more difficult to evaluate numerically. The integration over the energy

for the non-equilibrium term has to be computed along the real axis, which is numerically

noisy due to the singularities and resonances arising from the three-dimensional 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 (4-20) can be equivalently

written as

I/'
', (r)= p dEp, (r) (4-23)
PlR

where

/PL
P~(r)= dE [p,(r)+p p(r)]+,p (r). (4-24)


In principle, the electron density evaluated from either equation (4-20) or equation

(4-24) 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' (4-25)


where









L 2
c neq,T (4-26)
SL R
[neq, +[pneq,a]

with

,L = f dEp-
Pneq,T d ,(-
PR (4-27)

= 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 [62-64]. 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 two-dimensional 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 self-consistent 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 self-consistent

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 4-1. 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. 4-1 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. 4-2 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 4-2. 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 two-dimensional (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. 4-2 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 4-3. 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.

4-2. 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. 4-3. 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. 4-3 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 (4-28)
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. 4-4, 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 (10-5)


.0 0.2

.5 0.15

S0.05


00'


Transmission (10-4)


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 4-4. 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 4-5. 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 4-6. 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 (4-22). The effect of the bias

voltage, through the change of the electro-chemical 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 non-self-

consistent calculations, and, conversely, how important is charge self-consistency 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 self-consistent 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 4-5 shows the I-V 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 (4-28) which is1.17 x 107 /(Qm2) Note also that the

non-self -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 self-consistent 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 4-7. 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. 4-5.

Figure 4-6 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 self-consistent zero-bias 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. 4-7 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 I-V curve for the

majority spin channel may extend up to 0.108 V. For the antiparallel case, Fig 4-8 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 4-8. 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. (4-29)
S R /'R










In Fig. 4-9 and 4-10, 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 4-7 and 4-8. 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 4-9. 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 (4-12) and (4-13) show the integrated transmission as a function of k// for

all k// points at a higher bias of 0.378 V. Compared to Fig 4-4, 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 linear-response limit is (Gp-

Ga)/Ga as given in the introduction. The voltage dependence of the TMR is shown in Fig.

4-13. 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 non-SCF calculation for TMR

is similar to that of the tunneling current for parallel spin alignment. This indicates that

the charge self-consistency is more important for the majority spin channel than the

























0 1 1---I 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 4-10. 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 (10-5)


0.5




Kx 0.1 0.


.0I


S0.2
0.15
~ 0.1
' 005
Ky


Figure 4-11. 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 (10-5)


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 4-12. 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. 4-13 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 voltage-dependence 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:Non-SCF' ----------

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 4-13. 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 spin-polarized calculation using the

non-equilibrium first-principle method for spin-dependent 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 non-SCF

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 self-consistent non-equilibrium calculation.














CHAPTER 5
PROBLEMS IN THE CURRENT NON-EQUILIBRIUM 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 first-principles 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 self-consistent procedure based on

DFT to determine the effective single-electron 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 single-particle transmission function over energy. Such

approaches have been used in many systems with varying degrees of success: the shape

of the I-V curve can be compared to the experiment and the structure of I-V 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 I-V 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 Hohenberg-Kohn 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

one-to-one 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 non-equilibrium 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 non-equilibrium open system. Then we propose a modification of the

Hohenberg-Kohn 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 non-equilibrium

open systems.









5.2 Hohenberg-Kohn 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. 1-1. 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. 1-1) 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

first-principles 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 Hohenberg-Kohn

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 Hohenberg-Kohn 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 Hohenberg-Kohn theorem must be

modified in order to be applied to such systems.

5.2.1 Voltage-dependence of the Density Functional

In order to see how the Hohenberg-Kohn theorem should be changed, we start from

a rigorous approach to this problem based on the adiabatic time-dependent procedure. We

take the same time-dependent 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 time-dependent system can be obtained by solving the

time-dependent Schrodinger equation

.0
S (t) = H(t)(t). (5-1)

The Hamiltonian can be written as

H(t) = H +V(t)

IH,= 1V,2+i + Vext(r,)+Z Vee(rj) (5-2)
2 1 I V(t)= V(t)( ect)


where, in the Hamiltonian, the first term H, consists of the kinetic energy, the nuclear-

electron interaction Vex and the electron-electron interaction V,. The second term is an

extra time-dependent 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 time-dependent system, Runge and Gross [69] showed that the wave

function at time t, ry(t), can be determined by the time-dependent electron density p(t)

and the initial state, up to a time-dependent phase factor,

,(t) = e 'O't)A/P(O, /(-o0)] (5-3)

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 Hohenberg-Kohn theorem for ground state, it is completely determined by the

initial electron density; hence ((-oc) can be eliminated in the equation (5-3), so that

y/(t) is determined up to a phase factor by p only.

Neglecting the phase factor, we come to the time-dependent DFT in this case,

which claims that any time-dependent physical quantity can be expressed as the

functional of time-dependent 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 Hohenberg-Kohn 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 (5-3) can be written as,

f(t) = e p(t)/[P(t),1 ] (5-4)

which means that at a given Vb, ry(t) is determined by p to a time-dependent 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

steady-state 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 exchange-correlation energy

functional Ex, [].

The voltage-dependence of those density functionals actually comes from the

voltage-dependence of the initial state. The first Hohenberg-Kohn theorem for the

equilibrium state is just a special case forVb = 0.

The voltage-dependence 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 voltage-dependence of the initial state

varies with the system, the voltage-dependence 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 Non-equilibrium Electron Density

The key to eliminate the voltage-dependence of the functional in the equation (5-4)

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 non-equilibrium system, the electron density can be expressed as [20, 22],


p(r)= I'd (r, E) (5-5)
21i

where, G'(r, ) is the non-equilibrium 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 non-equilibrium part p,,

p(r)) (r)+ pn(r) (5-6)

where,

1 ^R
2ii
P,, (r)=- 2IP dsG(r,E)
7 (5-7)

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) (5-8)
-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, ). (5-9)
i"
PLR

Since the left lead is in local equilibrium, the first Hohenberg-Kohn 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 (5-8),


we see that PRt changes to PR + Vy. Introduce the new pR and GR (E) into the equation

(5-9), 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 (5-4) can be rewritten as,

Y(t =0) = e '(t =o[pq (t = 0), p (t = 0)] (5-10)









where, the voltage-dependence 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 time-dependent system, we can define the action A as,


A = ~(t)i- IH(t) V(t))dt (5-11)


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 (5-12)
9(0

is completely equivalent to equation (5-3), 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 (1-1). Thus, the Equation (5-12) can not be directly applied.

However, for the case of time-dependent 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 (5-12) to calculate the action and the stationary requirement becomes,

9A
= 0 (5-13)
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 non-equilibrium









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 (5-14)
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 time-dependent process, the action reduces to

0
A= E(t)dt. (5-15)


Then, the Equation (5-14) becomes,

= 5E0 (5-16)
3p,q (t) gp, (t)

where, E is the total energy of the system at time t.

Both the ground state and the time-dependent DFT work in the framework of a

single-particle picture. We follow the same procedure. Consider an independent particle

system, which has a same geometry as shown in Fig. 1-1 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 well-separated 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 time-dependent 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 non-interacting, both Green's functions G (E) and G' () are

diagonalized by m,~ It is straightforward that equation (5-7) becomes,


Peq= Z ,e 2
m2 (5-17)
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 (5-18)


Assume that Peq and p,, yielded by equaton (5-17) and (5-18) 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 single-particle kinetic

energy,

T = T +T

Teq Z (O -V2 m). (5-19)
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 (5-20)
2-: j- 2









where EcoU, is the classical electron-electron 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 .(5-21)

Ve" = E
ef Pn

The Equation (5-16) can be written in terms of those orbitals as,


(E V2 -), R = 0
2 (5-22)
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 (5-22) 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 non-equilibrium electrons see a different mean field potential than the equilibrium

electrons do. The self-consistent procedure makes two electron densities Peq and p, to

converge respectively, and the current is carried by the non-equilibrium electrons.

5.3 Exchange energy for uniform electron gas

The exchange-correlation 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 single-particle orbitals can be written as,


(k, k, k) = e (kx+kk ) (5-23)
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

right-going states (k > 0) and the left-going states (k < 0). From equations (5-17) and

(5-18), we see that the occupation of those states is: belowke, all states are occupied and

between k, and k,, only right-going 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 right-going 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 (5-24)
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 (5-25)
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 non-equilibrium 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 (5-26)
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. 5-1 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 voltage-dependence 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 5-1. Exchange energy at different bias voltages for different p (Electrons/m3)

We notice that although the voltage-dependence 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. (5-26) as,









K =Ke +K +K

Keq f Ieq

1 :1 2 (5-27)
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 non-equilibrium electrons and the third one denotes the

contribution of exchanging between equilibrium and non-equilibrium electrons. The

voltage-dependence of these three terms for p = 7.29 Electrons/m3 is shown in Fig. 5-2.

At bias zero, since p, is zero, there are no non-equilibrium 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 (5-28)

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. 5-3.








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 non-equilibrium 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 5-2. Split of the exchange energy: (a) the exchange between non-equilibrium






-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


(5-29)


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 5-3. 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


(5-30)


where, C, = 0.7386.


Equation (5-30) gives the exchange energy functional of the uniform system for all

bias voltages. It clearly shows the dependence on the non-equilibrium electron density as

well as the dependence on the total electron density, which support the analysis of the


-m-p=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 non-equilibrium

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 I-V curve can be compared to the

experiments and peaks in the curve can be explained by tunneling through molecular

orbitals.

We apply the state-of-art approach which combines the non-equilibrium 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 light-driven molecular switch,

which is important for application of nanoelectronics. The I-V 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 bias-enhanced 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 non-equilibrium 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 voltage-dependence of exchange energy.

By introducing 'non-equilibrium electron density', we show that the voltage-dependence

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 DE-FG02-02ER45995.














APPENDIX A
EXPANSION OF LESSER GREEN'S FUNCTION

The non-equilibrium Green's functions are defined as

G(1,1') = -i( )(1 (A-l)

where, the shorthand notation 1 (x, t,) is used. The time-order 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 non-equilibrium Green's

functions:

G(1,) t, t, e+
G'(1,1') tl e-,t, e+
G(J1,') = (A-2)
G'(1,1') tj +, t,e-
G' (1,1') ti, ti, -

These four Green's functions can be denoted as the time-ordered 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 (2-1), a particular lesser Green's

function Gk defined as


G (, (A-3)






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 A-1. 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,,,, (A-4)


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 (A-4) can be used for both equilibrium and non-

equilibrium Green's functions. For non-equilibrium case, since the time of the

perturbation can be either on + branch or branch, the equation (A-4) gives two terms.







87


As an example, the expansion of the lesser Green's function Gk contains two terms as


shown in Fig. A-1. Thus, Gk can be written as


Gkn = 7Vkam kakaGn gakaGmn. (A-5)
m














APPENDIX B
DETAILS OF NON-EQUILIBRIUM SELF-CONSISTENT 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 self-energies of two leads according

to equation (3-2). Since the tight-binding model is used to describe two leads, only the

part L (R) of the device region is affected by the left (right) semi-infinite lead. Thus the

dimension of the self-energy 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, self-energies of two leads are not affected by the device region, so they

don't change in the self-consistent procedure.

The self-consistent 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), (B-l)

2- P-r

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', = (B-2)
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 (B-2) can be expanded as

K+SLL HLL HL +SML HML E+SRL HRL
G' = E '+S HLM ES -Ha ~+ Sm H -H (B-3)
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 (2-11)

and the density matrix is determined by


jP, = IG; (8)dE, (B-4)

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(), (B-5)
1,j 1,j

where, 0, (i) is the Gaussian basis in real space.

The electric static potential is defined as

S() () + Vb, () (B-6)


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') (B-7)
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.




Full Text

PAGE 1

TRANSPORT PROPERTIES AT NANO SCA LES VIA FIRST PRINCIPLES STUDIES By CHUN ZHANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004

PAGE 2

Copyright 2004 by Chun Zhang

PAGE 3

To my parents.

PAGE 4

ACKNOWLEDGMENTS First, I would like to thank my advisor, Dr. Hai-Ping 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 long-term 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. iv

PAGE 5

TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................iv LIST OF FIGURES ..........................................................................................................vii ABSTRACT .........................................................................................................................x CHAPTER 1 OVERVIEW.................................................................................................................1 2 FIRST PRINCIPLES APPROACH COMBINING NGFT AND DFT........................4 2.1 Model Hamiltonian and Landauer-type Formula...................................................4 2.2 The Combination of DFT and NGFT for An Open System...................................8 2.2.1 Self-Consistent Procedure for Ground State of Closed System...................8 2.2.2 Self-consistent Procedure for Open System...............................................10 3 COHERENT ELECTRON TRANSPORT THROUGH AN AZOBENZENE MOLECULE: A LIGHT-DRIVEN MOLECULAR SWITCH..................................12 3.1 Introduction...........................................................................................................12 3.2 Tight-binding Model and Self-energies................................................................14 3.3 Conductance at Zero Bias.....................................................................................18 3.4 The Switch under A Finite Bias............................................................................23 4 QUANTUM TUNNELING THROUGH MAGNETIC JUNCTIONS......................30 4.1 Introduction...........................................................................................................30 4.2 Non-equilibrium Greens Function Formalism in Real Space.............................34 4.3 Self-consistency 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 NON-EQUILIBRIUM APPROACH FOR OPEN SYSTEMS BASED ON DFT.....................................................................................63 5.1 Introduction...........................................................................................................63 5.2 Hohenberg-Kohn Theorem for Open System under Finite Bias..........................65 v

PAGE 6

5.3 Exchange Energy for Uniform Electron Gas........................................................75 5.4 Conclusion............................................................................................................82 6 SUMMARY AND CONCLUSIONS.........................................................................83 APPENDIX A EXPANSION OF LESSER GREENS FUNCTION.................................................85 B DETAILS OF NON-EQUILIBRIUM SELF-CONSISTENT PROCEDURE...........88 LIST OF REFERENCES...................................................................................................90 BIOGRAPHICAL SKETCH.............................................................................................94 vi

PAGE 7

LIST OF FIGURES Figure page 1-1. Partition of the open system connected to two metal electrodes..................................2 3-1. Physical models of an azobenzene switch in contact with two gold leads: trans isomer (upper) and cis isomer (lower).....................................................................14 3-2. Tight-binding model of an open system connected to two semi-infinite leads. Indices denote nodes of device region and are the nearest node of left and right lead. Only the nearest-neighbor interaction is considered.......................................15 ,ij ,ab 3-3. A perfect lead that consists of two semi-infinite leads and a hopping term. Indices m and i belong to the left semi-infinite lead and i+1, j belong to the right semi-infinite lead...........................................................................................................................16 3-4. 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..............................................................................................19 3-5. 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 3-6. Local DOS at the Fermi energy projected on seven sites in the switch. From 1-7, or from left to right, Au lead, S atom, benzene ring, double-bonded N 2 another benzene ring, S atom and Au lead............................................................................22 3-7. Non-equilibrium self-consistent procedure. The non-equilibrium Greens function is calculated according to equation (2-11)...................................................................24 3-8. The difference of electronstatic potential between bias one volt and bias zero for cis configuration. The x-z plane is parallel to the N2 bond, cutting through the midpoint of Au unit cell...........................................................................................25 3-9. The transmission probability as a function of single-particle energy and the external bias voltage for cis configuration.............................................................................26 3-10. 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 vii

PAGE 8

3-11. The differential conductance as a function of bias voltage for both trans and cis configurations...........................................................................................................28 3-12. The transmission as a function of energy at 0.64, 0.72 and 0.8 Volts for cis configuration............................................................................................................28 4-1. Shifts in the density of states in each atomic layer at a bias voltage of 0.54 V..........44 4-2. 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 4-3. 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..............................................................................47 4-4. 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 4-5. 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 4-6. 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 4-7. Transmission probability at energy RE as a function of k // along the line of for the majority spin channel for parallel alignment of the moments at different bias voltages..............................................................................................53 xkk y 4-8. Transmission probability at energy RE as a function of k // along the line of for the minority spin channel for antiparallel alignment of the moments at different bias voltages..............................................................................................55 xkk yyy 4-9. Integrated transmission probability as a function of k // along the line of for the majority spin channel for parallel alignment of the moments at different bias voltages.....................................................................................................................56 xkk 4-10. Integrated transmission probability as a function of k // along the line of xkk for the majority spin channel for parallel alignment of the moments at different bias voltages.....................................................................................................................58 viii

PAGE 9

4-11. 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 4-12. 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 4-13. 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 5-1. Exchange energy at different bias voltages for different (Electrons/m 3 )..............78 5-2. Split of the exchange energy: (a) the exchange between non-equilibrium electrons and the exchange between non-equilibrium and equilibrium electrons and (b) the exchange between equilibrium electrons.................................................................80 5-3. Exchange energy as the functional of /n for different electron densities............81 A-1. Expansion of lesser Greens function. The time of the perturbation can be either on + branch or on branch............................................................................................86 ix

PAGE 10

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: Hai-Ping 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 many-body 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 non-equilibrium 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 Non-equilibrium Greens 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 many-body effects and NGFT naturally includes the non-equilibrium effects. The new approach uses NGFT to treat the non-periodic 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 light-driven molecular switch. The switch, made of a single molecule, is one of the most x

PAGE 11

important elements of nano-electronics. However, most proposed molecular switches are driven either by an external bias voltage or by STM manipulation, neither of which is ideal for nano-scale circuits. The switch we designed has a high on-off conductance ratio and more importantly, can be driven by photons. In following chapters, we generalize the method to the spin-dependent case and apply it to a magnetic layered structure. We implemented the method within the framework of the Layer Korringa-Kohn-Rostoker (LKKR) approach, which is particularly well-adapted to the layered structure and found a bias-enhanced tunneling magneto-resistance (TMR) for the Fe/FeO/MgO/Fe junction. Our results are important not only for application, but also for understanding of the voltage-dependence 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 impurity-free layered structure, a different behavior of TMR may occur. xi

PAGE 12

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 modern technologies, such as STM and AFM, to make molecular transistors and switches etc. [1-5]. 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 time-dependent 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 steady-state 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 steady-state 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 1

PAGE 13

2 conjunction with the effective-mass assumption. The Boltzman equation is based on a probabilistic description for the evolution of a single-particle distribution which is a function of classical variables ),,(tprf r and p However, when the device region is at nano-scale, 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 non-equilibrium Greens functions technique (NGFT) and Density Functional Theory (DFT) has proven to be successful [6-22]. The method can be summarized in two steps. First, a non-equilibrium self-consistent procedure based on DFT and NGFT is used to determine the single-particle effective potential. Second, a Landauer-type formula is used to calculate the current, which relates the current to the integral of the single-particle transmission probability over energies. Current I Left lead Right lead Device Region Figure 1-1. Partition of the open system connected to two metal electrodes. The approach combining NGFT and DFT is based on the single-particle 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

PAGE 14

3 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 two-step 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 light-driven molecular switch. In Chapter 4, we generalize the method to spin-dependent tunneling junctions with application to a layered magnetic Fe/FeO/MgO/Fe junction. The calculation suggests a bias-enhanced tunneling magneto-resistance (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 two-step approach. For the first step, in which the effective single-particle potential is determined, DFT is the biggest problem. The Hohenberg-Kohn 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 non-equilibrium effects into DFT. For the second step, in which a Landauer-type formula is used to calculate the current, the problem is that the Kohn-Sham orbitals are used to calculate the current. Such orbitals are not genuine physical quantities, which causes trouble in choosing basis sets and functionals.

PAGE 15

CHAPTER 2 FIRST PRINCIPLES APPROACH COMBINING NGFT AND DFT 2.1 Model Hamiltonian and Landauer-type 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 non-interacting, 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 CouCenRLHHHHH (2-1) where LH RH and CenH denote the Hamiltonian for the left lead, the right lead, and the device region respectively. The last term CouH denotes the coupling between leads and the device region. Since the two leads are non-interacting, we have nRLknknkCounnnCenRLkkkkRLchdcVHddHHccH,,,int,,,,,,.].[}){},({ (2-2) 4

PAGE 16

5 where indices ,k denote the channel of left and right leads (Bloch wave vectors and spin) and denotes the single-particle basis of the device region. The key issue is how to treat the interaction of the device region. n The expression for the current can be derived by taking the time derivative of the total number of electrons in the left lead LLLNHieNeJ, (2-3) where and is the current from the left lead into the device region. Inserting the Hamiltonian into the expression for the current, we have LkkkLccN, LJ nLknknkknnknLkknnknknkLGVGViecdVdcVieJ,,*,,,,*,, (2-4) where the current is expressed in terms of two non-equilibrium lesser Greens functions. Now, the problem becomes how to calculate the lesser Greens functions. This can be achieved by a perturbation procedure. At t 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 R is higher than the Fermi energy of left lead L As we gradually turn on the coupling from to t 0 t the system will achieve a steady state at 0 t with a current from left to right. Taking the coupling term in the Hamiltonian to be the perturbation, the lesser Greens functions can be expanded as J

PAGE 17

6 mamnkkmnrkkmknkmakknmkkrnmmkknGgGgVGgGgGVG][][,,,,,,*,, (2-5) where runs over all basis states in the device region ,denotes the Greens function for non-interacting leads and the upper indices m g a r denote retarded and advanced Greens functions respectively. The details of the expansion can be found in Appendix A. Introducing the equation (2-5) into equation (2-4) and doing a Fourier transform gives ]})()()()([Re{22,,,*,, nmakkrnmkkmnLkmknkLGgGgVVdeJ (2-6) where the summation runs over all channels of the left lead. Since the lead is non-interacting, the Greens functions can be calculated as agg, igifgkakkkkkk1)()()(2)(,, (2-7) where )( kf is the Fermi-Dirac distribution function. Insertion of the equation (2-7) into equation (2-6) gives the current as ]})([{2GGGfTrdieJarLLL (2-8) where the line-width function L is defined as LkmknkmnkLVV)()()(2)]([* (2-9) in which the term is the density of states in the lead. A similar formula can be derived for, the current from the right lead into the device region, so the total current is given by RJ )(21RLJJJ which is

PAGE 18

7 ]}][[]{[22arRRLLRLGGffGTrdieJ (2-10) where the current is related to the local properties of the device region. In equation (2-10), the Greens functions in the device region are solvable for only a few models. The simplest one among them is the freelectron model, which means that the device region is also non-interacting. For the free-electron model, we have [15, 20, 25] ][RRLLarffiGGG (2-11) and then the current can be written as }{][LrRaRLGGTrffdeJ (2-12) which is a Landauer-type formula. Since the interaction in both leads and the device region usually plays an important role in the tunneling process at nano scales, equation (2-12), 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 single-particle picture. Within present-day approximations, it is able to take into account most many-body effects. Moreover, due to the success of DFT in the electronic structure calculation of molecules, it is natural to combine DFT and the Landauer-type formula to study molecular electronics.

PAGE 19

8 2.2 The Combination of DFT and NGFT for An Open System 2.2.1 Self-Consistent 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 [26-30]. The main ideas of DFT are outlined in the following paragraphs. The basis of DFT is the Hohenberg-Kohn 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 rdrVrFEext )()(][][ (2-13) where the term represents the entire external potential (including the ion-electron interaction) and the first term extV ][ F is a universal functional related to the electron-electron interaction and the kinetic energy of electrons. The theorem also shows that the true ground-state electron density minimizes the energy functional ][ E The only thing used in the proof of Hohenberg-Kohn Theorem is the minimum energy principle of ground state. On the basis of the Hohenberg-Kohn theorem, a self-consistent procedure based on the single-particle picture can be constructed. Assume that there is a non-interacting system which yields the same electron density as the interacting system of interest, then the electron density can be written in terms of the single-particle functions iirr2)()( (2-14)

PAGE 20

9 which are usually called Kohn-Sham orbitals. The kinetic energy of the non-interacting system can be calculated in terms of the Kohn-Sham orbitals (Hartree atomic units) iiisT2][2 (2-15) Then the exchange-correlation energy can be defined as ][][][][ HsxcETFE (2-16) where the term ][ HE is the classical Coulomb energy. The exact form of ][ xcE is unknown. However many kinds of exchange-correlation functionals, most based on the local density approximation, have been proposed and successfully used. By minimizing the energy functional, we get the Kohn-Sham equations )()()]([)()]([212rrrVrVrViiixcextH (2-15) where the first term is the well-known classical Coulomb potential and the term is the exchange-correlation potential defined by HV xcV ][xcE The effective Kohn-Sham potential is defined as )]([)()]([)]([rVrVrVrVxcextHeff (2-16) Now, a self-consistent procedure (Kohn-Sham procedure) can be completed as follows. Given the electron density the effective Kohn-Sham potential ][ effV can be calculated according to equation (2-16). Knowing ][ effV we can get the Kohn-Sham orbitals i by solving equation (2-15) and then the new electron density can be calculated according to equation (2-14). Repeat the procedure until the self-consistency is achieved. The Hohenberg-Kohn theorem guarantees that the Kohn-Sham procedure yields the true

PAGE 21

10 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 self-consistent procedure based on the single-particle picture. 2.2.2 Self-consistent Procedure for Open System For an open system with the geometry shown in Fig. 1-1, 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 Kohn-Sham 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 Kohn-Sham 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 self-consistently in the device region with two fixed, open boundaries. The Greens functions technique described below is an efficient way to do this. The electron density in the device region is directly related to the lesser Greens function by dGi)(21 (2-17)

PAGE 22

11 For a single-particle Hamiltonian, the lesser Greens function can be calculated by the equation (2-11), where the retarded and advanced Greens functions are related to the effective potential of the device region and two open boundaries as )()()(1arRarLCenarHEG (2-18) In equation (2-18), the term CenH denotes the effective single-particle Hamiltonian and the two terms are self-energies 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. RL, As a conclusion, the self-consistent 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 in the device region, the effective Kohn-Sham potential can be calculated. Knowing the effective potential, the lesser Greens function can be computed and then a new electron density is given by equation (2-17). After self-consistency is achieved, the Landauer-type formula (2-12) can be used to calculate the current in the system. The non-equilibrium self-consistent procedure mentioned above is implemented in a Greens functions representation which employs a multi-band tight-binding model in Chapter 3 and is implemented in a wave function representation which includes spin in Chapter 4.

PAGE 23

CHAPTER 3 COHERENT ELECTRON TRANSPORT THROUGH AN AZOBENZENE MOLECULE: A LIGHT-DRIVEN MOLECULAR SWITCH 3.1 Introduction Electronic devices soon will be constructed of molecules or have features of molecular size. The results of recent experimental [1-5] and theoretical studies [6-20] 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 [31-35]. 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, light-driven processes, are therefore highly desirable. The photosensitive azobenzene molecule discussed here is one such candidate for an opto-electronic device. The azobenzene molecule has two stable conformations [36-38], 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 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]. 12

PAGE 24

13 Previous work has shown that a nanoscale opto-mechanical 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, Au-S-azobenzene-S-Au. It consists of a single azobenzene molecule, which has been functionalized by replacing a hydrogen atom with a CH 2 S group to provide a good contact with gold, and two semi-infinite gold leads. Note that the name of the modified molecule is technically bis-(4-methanethiol-pheyl) 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. 3-1 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 and, for the cis configuration, 16.23 In subsequent sections, we study the transport

PAGE 25

14 properties for both equilibrium and non-equilibrium case. In both cases, the localized basis set and the tight-binding model is employed. Figure 3-1. Physical models of an azobenzene switch in contact with two gold leads: trans isomer (upper) and cis isomer (lower). 3.2 Tight-binding Model and Self-energies As mentioned in Chapter 2, the key issue is to evaluate the self energies of two semi-inifinite leads, and. Once these two terms are computed, we can construct the retarded and advanced Greens functions and of the device region in the presence of two leads. Then the lesser Greens function can be calculated according to equation (2-11). L R rG aG The self-energy term represents the effects of the semi-infinite 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 tight-binding model is a natural choice. Fig. 3-2 sketches the model, where indices j i, denote nodes of the device region and are the nearest node of the left and right lead respectively. Note that hopping terms are operators [42]. ba, ',hh

PAGE 26

15 i j abh'h Figure 3-2. Tight-binding model of an open system connected to two semi-infinite leads. Indices j i, denote nodes of device region and are the nearest node of left and right lead. Only the nearest-neighbor interaction is considered. ba, Within the framework of tight binding, the Greens function can be written as ijG jjbbCenijijaaCeniiCenijijGhghGhGghGGG'' (3-1) where the two hopping terms are taken as perturbations and the Dyson equation is used. The term represents the Greens function of the device region without two leads and the term denotes the unperturbed Greens function of the leads. Then, self-energies can be defined as CenG g ''hghhghbbRaaL (3-2) where, are self-energies of the left and right lead respectively. Two hopping terms are projections of the effective Kohn-Sham Hamiltonian to a localized basis of each node (hopping terms are matrices also). The only problem remaining is two so-called surface Greens functions and. RL, aag bbg Since a perfect lead with periodic boundary conditions consists of two semi-infinite leads and a hoping term within the framework of tight-binding as shown in Fig. 3-3, it is possible to construct the Greens functions of semi-infinite leads on the basis of Greens functions of a perfect lead by a perturbation procedure. The whole procedure is described in the following.

PAGE 27

16 h m i i+1 j Figure 3-3. A perfect lead that consists of two semi-infinite leads and a hopping term. Indices m and i belong to the left semi-infinite lead and i+1, j belong to the right semi-infinite 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 Greens function of the perfect lead and G 0 the Greens function of two semi-infinite leads. The Dysons equation reads qjpqipijijGVGjVGGiGG000 (3-3) where q p run over the whole system. Since V only contains two terms V i,i+1 and V i+1,i P equals to i or i+1. Because G 0 i,i+1 =0 (without perturbation, an electron can not propagate from i to i+1), we have jiiijiiiiiijhGgGVGG,1,11,0 (3-4) Similarly, (3-5) miiimiiiiimiGhgGVGG,1,1,,101,1,1 In order to solve g from equation (3-4) and (3-5), we have to know relations between G ij and G i+1,j This can be found by solving the Greens functions of the perfect lead. For a perfect lead, the tight-binding Hamiltonian is mmmmmmmmaahahaaahH110 (3-6) where are matrices, being the number of orbitals within one unit cell of the lead. The eigenvalue problem for Bloch states takes the form hh,0 NN N

PAGE 28

17 0),()()(*0kuehehhikdikd (3-7) where is the Bloch vector and is the lattice constant. The energy term k d contains an infinitesimal imaginary part, with iE 0 Equation (3-7) can be written as 02iiiuhhM (3-8) where, and It is to be noted that not all eigenvalues 0hM hh i of equation (3-8) for a given energy represent the Bloch state. Because equation (3-8) is quadratic in 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 i ) and another half are inside the unitary circle (we denote them as i ). Using the superscripts to denote eigenvalues and eigenvectors for and, respectively, we can introduce two matrices, NN and 1*1)()(UUUU (3-9) where is a matrix, whose th column is the eigenvector and is a diagonal matrix, whose th diagonal element is the eigenvalue The matrices, U i iu i i and satisfy the matrix equation 0)()(0)()(22hhMhhM (3-10) The retarded Greens function is defined by the matrix equation. We have IGH)( .0)()(,1,1jijiijijGhGhGMjGHi (for j i ). (3-11)

PAGE 29

18 For a perfect lead, the solution of this equation is done via recursion relations of Greens functions, and. This gives the equation jiijAGG,1 jijiGAG,12,1 0. (3-12) )(2AhhAM Comparing with equation (3-10), we see that, for j i Introducing into equation (3-5), we have A jiijAGG,1 1)(hgii (3-13) Similarly, for j i we have (3-14) 11,1)(hgii 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 linear-response 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], )(22FEThe (3-9) where 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, STO-3G Gaussian basis [28] for all other atoms (C, S, N, H). The pbe0 exchange-correlation functionals [26] are used. For the wire system depicted in Fig. 3-1, the system is, in principle, infinite, and periodic boundary conditions are not appropriate. However, due to )(FET

PAGE 30

19 screening, the Kohn-Sham potentials in the two semi-infinite 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]. -0.24-0.22-0.20-0.180.00.40.8 CisCisTransTransmission Density of states (in device area)Energy (Hartree)EfTrans0200400 aa'bb' Figure 3-4. 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 Greens function techniques, the effect of the two leads on the switch can be absorbed into self-energy terms, which are calculated using the same method as mentioned in previous section. The transmission function T is evaluated by the Caroli formula [25] )(E F ][aRrLGGTrT (3-10) which is the kernel of the integral in equation (2-12). The density of states (DOS) in the switch region can be calculated from the trace of the retarded Greens function, which

PAGE 31

20 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. -0.24-0.22-0.20-0.180200 Local density of statesEnergy (Hartree)Au layersaa'b'0200 Trans Sulpher atoms a0200 CisCisTransCis EfTransAzobenzene moleculebb' Figure 3-5. 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. 3-4. Unlike the DOS of an isolated molecule, Fig. 3-4 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

PAGE 32

21 and b, for the trans isomer, and peaks and for the cis isomer. These peaks are indications of localized states, which can be seen more clearly in the local density of states. 'a b As depicted in Fig. 3-5, 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 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 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. 3-5 indicates clearly that the broad background in the DOS derives primarily from the Au layers. a 'a Figure 3-5 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 zero-bias 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. 3-4. For the trans isomer, tunneling through the perturbed HOMO leads to a very broad and

PAGE 33

22 intense peak in the transmission function, but a sharp and smaller peak for the cis isomer (Fig.3-4, 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 zero-bias conductance for the trans structure is which indicates a moderately good conductor. In contrast, the conductance for the cis structure is nearly 2 orders of magnitude lower. 141021.0 1234567051015202530 Local density of statesAtom groups1: Left Au layer2: Left S atom3: Left benzene ring4: N=N5: Right benzene ring6: Right S atom7: Right Au layerCisTrans Figure 3-6. Local DOS at the Fermi energy projected on seven sites in the switch. From 1-7, or from left to right, Au lead, S atom, benzene ring, double-bonded N, another benzene ring, S atom and Au lead. 2 To understand the tunneling process at the Fermi energy, Fig. 3-6 depicts the local DOS at the Fermi energy, Fig. 3-6 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 CH 2 group, the double-bonded N 2 the right benzene ring with a CH 2 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

PAGE 34

23 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 non-equilibrium self-consistent procedure discussed in Chapter 2 to

PAGE 35

24 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 zero-bias case. We neglect the effects of the so-called current-induced force by making this assumption. dGi)(21 ][ VVeff Figure 3-7. Non-equilibrium self-consistent procedure. The non-equilibrium Greens function is calculated according to equation (2-11). The self-consistent procedure is depicted in Fig. 3-7. Since we are working within the framework of DFT, the most important physical quantity is the electron density. Given the electron density the effective potential can be computed as a functional of the density. Once we know the effective potential, the non-equilibrium Greens function can be calculated according to equation (2-11), where the retarded and advanced Greens functions are calculated in the same way as in the previous section. The details of the self-consistent procedure can be found in Appendix B. effV G One of the major accomplishments of the non-equilibrium first principles calculation is that the electrostatic potential due to the external bias can be computed self-consistently. Fig. 3-8 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 x-z plane is parallel to N double bond of cis configuration, cutting

PAGE 36

25 through the midpoint of the Au unit cell. The drop of the electrostatic potential across the switch area can be clearly seen. Figure 3-8. The difference of electronstatic potential between bias one volt and bias zero for cis configuration. The x-z plane is parallel to the N2 bond, cutting through the midpoint of Au unit cell. The current through the system is calculated by equation (2-12), where the current is expressed as the integral of the single-particle 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 3-9, we show the transmission probability as a function of both energy and bias voltage, where the transmission is expressed as ),(),(),(),(),(VEGVEVEGVETrVETaRrL (3-13)

PAGE 37

26 Figure 3-9. The transmission probability as a function of single-particle 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 Kohn-Sham 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. I-V curves through the trans and cis configuration are shown in Fig. 3-10. The differential conductance, defined as the derivative of the current with respect to the bias voltage, is shown in Fig. 3-11. 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

PAGE 38

27 than the cis for all bias voltages, which shows the high stability of the switch under a finite bias. Figure 3-10. 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. 3-10. 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. 3-11. This can be explained by consideration of Fig. 3-12, 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 (3-1), 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. 3-10 clearly shows such an example, where as the bias voltage

PAGE 39

28 Figure 3-11. The differential conductance as a function of bias voltage for both trans and cis configurations. Figure 3-12. The transmission as a function of energy at 0.64, 0.72 and 0.8 Volts for cis configuration.

PAGE 40

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. 3-11.

PAGE 41

CHAPTER 4 QUANTUM TUNNELING THROUGH MAGNETIC JUNCTIONS 4.1 Introduction Quantum transport in open systems far from equilibrium is an important and challenging problem [1-5] 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 [6-20] with various degrees of success. The preceding chapter is an example. The problem of spin-dependent tunneling in magnetic heterostructures, however, presents a significantly greater challenge than the spin-independent case to the theory of non-equilibrium transport because one must simultaneously deal with the complexities of non-equilibrium transport and complex electronic structures along with magnetism. Spin-dependent tunneling junctions are considered to be the most likely candidate for next-generation devices for spintronics. Compared to the trilayer or multilayer giant magnetoresistance (GMR) systems, magnetic tunneling junctions tend to have a relatively higher magnetoresistance ratio [45-48], which is an important factor in device designs. The magnetoresistance ratio is defined as (G p -G a )/G a where G p is the conductance for the case that the alignments of net magnetic moments of two leads are parallel and G a is the conductance for the case that the alignments of net magnetic moments of two leads are anti-parallel. One of the 30

PAGE 42

31 significant differences between a GMR system and a spin-dependent 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 first-principles studies [50-54] 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 linear-response 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 [51-52]. 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 voltage-dependence of TMR in terms of the impurity scattering at the interface. They suggested that with a high quality barrier, a different voltage-dependence of TMR may occur. Some groups [55-57] 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

PAGE 43

32 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, first-principles calculation [6-20] of electronic structure and tunneling conductance under a finite bias voltage has been successfully applied to molecular electronics. As mentioned in Chapter 1, the first-principles approach can be built on either non-equilibrium Greens 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 Korringa-Kohn-Rostoker (KKR) approach [61], which was used successfully in calculating linear-response transport properties of spin-dependent 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 two-terminal conductance within the linear-response 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 [62-65]. The inclusion of a single atomic layer of FeO on one of the interfaces is important. An earlier

PAGE 44

33 linear-response 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 first-principles 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 exchange-correlation 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 exchange-correlation potential. In the following section, we will show that our approach is equivalent to the Keldysh Greens function approach within the single-particle picture, but without the

PAGE 45

34 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 Greens function. 4.2 Non-equilibrium Greens Function Formalism in Real Space The Keldysh formalism incorporates the open boundary conditions of a non-equilibrium system into the components of the Greens functions. Because not all of the Keldysh Greens functions are independent of each other, we can use the usual advanced and retarded Greens functions along with one additional component of the Keldysh Greens functions. Common practice is to use the lesser Greens function which is defined as G ),(),(),;,(''''trtritrtrG (4-1) where the angle brackets indicate the ensemble average and are field creation and destruction operators of electrons. There is a simple relation between the lesser Greens function and the electron density ),;,(),(trtriGtr Since we are only interested in the steady state solution, the time dependence is only through. So, the lesser Greens function can be Fourier transformed into 'ttt tdetrrGErrGdEeErrGtrrGtiEtiE);',(21);',();',();',( (4-2) Note that writing the Greens function with a single energy argument implies that the scattering in the tunneling region is elastic. The charge density is also independent of time, dEErrGir),,(21)( (4-3)

PAGE 46

35 The problem remains to evaluate );',(ErrG For non-interacting 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 Greens function in the form tdetrttriErrGtiE),(),(2),',(' (4-4) In the framework of mean-field theory, these two operators can be expanded as single-electron operators nnnnnnCtrtrCtrtr),(),(),(),(* (4-5) where n is the stationary state single particle wave function and two terms are creation and destruction operators for this stationary state. nnCC, Substituting equation (4-5) into equation (4-4), we have nnnnnnnnnnnrrfEEirrCCEEiErrG)()()()()()(),,('*'*' (4-6) where the product is the number of particles that occupy the state nnCC n and its ensemble average is the distribution function Equation (4-6) is a general expression for the lesser Greens function in terms of the stationary wave functions and the nf

PAGE 47

36 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 1-1. 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 )(LEf which is diagonalized by scattering function incident from the left lead, and likewise for the right lead distribution function )(REf and the scattering wave functions from the right lead. Both and are Fermi-Dirac distribution functions defined with their respective chemical potentials f f L and R Next, using the subscript l to denote all degenerate indices of and correspondingly for, equation (4-6) can be written as (neglecting the bound states for the moment) m lmmmmlllrrEEifrrEEifErrG)()'()()()()(),',(*'* (4-7) We compare this result with that obtained in Chapter 2. Rewrite equation (2-11) as aRraLrGGfGGfG (4-8) Both equations are valid within the single particle picture. However, equation (4-7) is more general than the equation (4-8), 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 (4-7) applies to the entire space whereas equation (4-8) can be applied only within the device region.

PAGE 48

37 Substituting equation (4-7) into equation (4-3), we obtain the scattering part of the electron density mmRmllLlscatEfEfr22)()()( (4-9) Including the contribution from the bound states, the lesser Greens function contains an additional term where and bbbbbrrfEEi)'()()(* bf b are the distribution function and the single-particle functions of bound states. Accordingly, the electron density also has an extra contribution from bbbbbEf2)( (4-10) The determination of is not unique and is still being investigated by researchers. bf 4.3 Self-consistency and Tunneling Current for A TMR Junction In this section we outline our approach for non-equilibrium self-consistent calculations of electronic structure of a TMR junction. Although the approach applies generally to any lead-device-lead ensemble, the particular Fe/FeO/MgO/Fe tunneling junction in our calculation has a layered structure with a two-dimensional periodicity along each layer. This approach allows a two-dimensional Fourier transform that yields a quasi one-dimensional problem for each reciprocal vector k // The left and right leads are semi-infinite in the direction perpendicular to the layers. They maintain different chemical potentials L and R Without loss of generality, we may assume RL 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.

PAGE 49

38 Within each electrode, the Bloch eigenstates can be divided into right-going states and left-going 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, 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 ,,,'''mmmlsllllltt device region (4-11) Left lea d Right lead and ,,,'''''lllmsmmmmmtt device region (4-12) Right lead Left lea d where are the transmission coefficients and are the reflection coefficients. The superscripts indicate the direction of travel along the z direction, and the indices label Bloch states in the left and right electrodes, respectively, and span over all possible values of k ',tt ',rr ml, // and the spin. If one takes the jellium approximation in both leads, the Bloch wave functions become plane waves and 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.

PAGE 50

39 At zero temperature, both distribution functions become the step functions )(,RLE Then from equations (4-9) and (4-10) we obtain )()()()(rrdErdErbEErl (4-13) where mmmElllEEErEEr22)()(,)()( (4-14) These equations hold for all regions including the device region and both leads. We can rewrite equation (4-13) as LRrdErrEReq)()()( (4-15) where the equilibrium part of the electron density corresponding to R is defined as RrrrdErbEEReq)()]()([)( (4-16) Note that the superscript R indicates the chemical potential R used in the calculation of yet the density is defined for all three regions. Req Req Now, let us consider the electron density in the two leads separately. In the right lead, the first term in equation (4-15) 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 within the right lead E

PAGE 51

40 lmmmlmErtEEr2)()()( (4-17) The electron density in the left lead contains more terms than the right lead. Again, neglecting the interference terms, we can write '2'''2)()()()()(lllllllllErREErEEr (4-18) The normalization of the Bloch wave functions is such that the integral over the volume of a unit cell llllllEEEErd223)()( equals to the density of states at energy E of the left lead, and likewise for the right lead. The self-consistent procedure for the non-equilibrium 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 semi-infinite (lead) parts are calculated using the bulk self-consistent 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].

PAGE 52

41 In order to see how the spin-dependent calculation is done, we write the electron density in equation (4-14) in terms of spin explicitly ////////////2,,,2,,,)()()()(kkkEkkkEEErEEr (4-19) where the term is the spin index. The subscripts now comprise both k ml, // and Remember that sum over ,E ll, and sum over ,E mm, Then equation (4-15) can be written as LRrdErEReq)()(,, (4-20) where the equilibrium part is )()()()(,,,rrrdErbEEReqR (4-21) The last term means only those bound states which have spin are counted. The equation (4-21) 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 spin-density functional to calculate the effective potential. Similarly, the current density for each spin is given by a Landauer-type formula lmlmlmtvvdEeJ,2, (4-22) where are the group velocities of the Bloch states, and only the Bloch states for the lmvv, 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.

PAGE 53

42 The equilibrium term in equation (4-20) can be evaluated using the standard contour integration technique for Greens function methods. The non-equilibrium term on the other hand, is more difficult to evaluate numerically. The integration over the energy for the non-equilibrium term has to be computed along the real axis, which is numerically noisy due to the singularities and resonances arising from the three-dimensional 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 Greens function directly, the integration still would need to be done along the real axis because the lesser Greens 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 (4-20) can be equivalently written as (4-23) LRrdErELeq)()(',, where )()()()(,,,rrrdErbEELeqL (4-24) In principle, the electron density evaluated from either equation (4-20) or equation (4-24) 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]: ')1(ccnew (4-25) where

PAGE 54

43 2,2,2,RneqLneqLneqc (4-26) with LRLRERneqELneqdEdE,,,, (4-27) 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 [62-64]. 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 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 while the distance between Fe and O atoms in the FeO layer is taken to be 2.1 Because we are unable to incorporate the considerations of chemical or structural disorder into the Keldysh Greens 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 two-dimensional 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 self-consistent part of the calculation, we use one energy point for biases below 0.01 Hartree and two energy points

PAGE 55

44 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 self-consistent 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 -0.2 -0.1 0 0.1 0.2 0.3 Shift in DOS (eV) Fe O MgO Fe Mg FeO Majority Minority Figure 4-1. 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

PAGE 56

45 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. 4-1 for parallel alignment of the moments. The DOS for the iron layers in the leads is shifted by roughly 2/1 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. 4-2 we show the capacitance 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 dVdQ/

PAGE 57

46 0 0.002 0.004 0.006 0.008 0.01 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 dQ/dV (electrons/volt) Bias Voltage (volts) Net charge Left electrode Right electrode Figure 4-2. 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 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 two-dimensional (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 Q

PAGE 58

47 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. 4-2 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 (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 9.22n 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 dm/dV ( /volt) B Bias Voltage (volts) Left Right Bulk electrode (left) Figure 4-3. 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.

PAGE 59

48 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. 4-2. 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 of Fe/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. 4-3. 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. 4-3 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

PAGE 60

49 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 lmmllmtvveG22 (4-28) 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. 4-4, 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.

PAGE 61

50 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Kx -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Ky 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Transmission (10 ) -5 (a) -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Kx -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Ky 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 Transmission (10 ) -4 (b) Figure 4-4. 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.

PAGE 62

51 -40 -30 -20 -10 0 10 20 30 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Current density (10 A/m ) 4 2 Bias Voltage (Volts) Majority Minority Majority:Non-SCF Figure 4-5. 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. -10 -8 -6 -4 -2 0 2 4 6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Current density (10 A/m ) 4 2 Bias Voltage (Volts) Current:SCF Current:Non-SCF Figure 4-6. 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.

PAGE 63

52 The tunneling current is calculated through equation (4-22). The effect of the bias voltage, through the change of the electro-chemical 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 non-self-consistent calculations, and, conversely, how important is charge self-consistency 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 self-consistent zero bias potential, in which the electrostatic potentials in the left and right electrodes are rigidly shifted by 2/1 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 4-5 shows the I-V 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 of. This number agrees very well with )/(1015.127m the result obtained from equation (4-28) which is. Note also that the non-self consistent calculation yields the same results within the linear region (for bias )/(1017.127m

PAGE 64

53 voltage up to about 0.1 V) but deviates significantly at higher biases, reaching about twice the self-consistent results at 0.5 V. The current of minority spin channels actually 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Transmission (10 ) -5 Kx (Ky) .0_Volts .054_Volts .108_Volts Figure 4-7. Transmission probability at energy RE as a function of k // along the line of for the majority spin channel for parallel alignment of the moments at different bias voltages. yxkk shows an oscillatory behavior. Since it is negligibly small compared to the majority spin channels, we cannot see this in Fig. 4-5. Figure 4-6 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 self-consistent zero-bias electrostatic potential. For most of the bias voltages the tunneling currents agree between the SCF and non

PAGE 65

54 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 RE and along a line of yxkk in the 2D reciprocal space, under different positive bias voltages. For the case of parallel alignment of the moments in both leads, Fig. 4-7 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 I-V curve for the majority spin channel may extend up to 0.108 V. For the antiparallel case, Fig 4-8 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

PAGE 66

55 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. 0 0.5 1.0 1.5 2.0 2.5 3.0 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Transmission (10 ) -5 Kx (Ky) .0_Volts .054_Volts .082_Volts .108_Volts Figure 4-8. Transmission probability at energy RE as a function of k // along the line of for the minority spin channel for antiparallel alignment of the moments at different bias voltages. yxkk 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 LRdEktkTRL)(1)(//// (4-29)

PAGE 67

56 In Fig. 4-9 and 4-10, we show this integrated transmission probability as a function of k // along a line of k x =k y in the 2D reciprocal space. The patterns shown in these two figures closely resemble those in Fig 4-7 and 4-8. 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. 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Integrated transmission (10 ) -5 Kx (Ky) .0_Volts .054_Volts .108_Volts Figure 4-9. Integrated transmission probability as a function of k // along the line of for the majority spin channel for parallel alignment of the moments at different bias voltages. yxkk 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

PAGE 68

57 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 (4-12) and (4-13) show the integrated transmission as a function of k // for all k // points at a higher bias of 0.378 V. Compared to Fig 4-4, 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 papJJJ/)( where are total currents for the parallel and antiparallel case respectively and its linear-response limit is (G apJJ, p -G a )/G a as given in the introduction. The voltage dependence of the TMR is shown in Fig. 4-13. 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 027.0 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 non-SCF calculation for TMR is similar to that of the tunneling current for parallel spin alignment. This indicates that the charge self-consistency is more important for the majority spin channel than the

PAGE 69

58 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Integrated transmission (10 ) -5 Kx (Ky) .0_Volts .027_Volts .054_Volts .082_Volts .108_Volts (a) 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Integrated transmission (10 ) -5 Kx (Ky) .108_Volts .136_Volts .163_Volts (b) Figure 4-10. Integrated transmission probability as a function of k // along the line of for the majority spin channel for parallel alignment of the moments at different bias voltages. yxkk

PAGE 70

59 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Kx -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Ky 0 0.5 1.0 1.5 2.0 2.5 3.0 Transmission (10 ) -5 Figure 4-11. Integrated transmission probability as a function of k // for the majority spin channel for parallel alignment of the moments at bias 0.378 V. -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Kx -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Ky 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Transmission (10 ) -5 Figure 4-12. Integrated transmission probability as a function of k // for the minority spin channel for antiparallel alignment of the moments at bias 0.378 V.

PAGE 71

60 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. 4-13 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 voltage-dependence 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

PAGE 72

61 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 TMR Bias Voltage (Volts) TMR:SCF TMR:Non-SCF A Figure 4-13. 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 spin-polarized calculation using the non-equilibrium first-principle method for spin-dependent tunneling junction under a

PAGE 73

62 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 non-SCF 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 self-consistent non-equilibrium calculation.

PAGE 74

CHAPTER 5 PROBLEMS IN THE CURRENT NON-EQUILIBRIUM 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 first-principles 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 self-consistent procedure based on DFT to determine the effective single-electron 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 single-particle transmission function over energy. Such approaches have been used in many systems with varying degrees of success: the shape of the I-V curve can be compared to the experiment and the structure of I-V 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 I-V 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 63

PAGE 75

64 molecular junction. Unfortunately, the problem is more fundamental than just the functional. The basis of DFT is the first Hohenberg-Kohn theorem [23], which states that for the ground state, the external potential is determined, within an additive constant by the electron density )(rVext 0V )(r Since all ground state properties are determined by the Hamiltonian, which is uniquely defined once the external potential is determined, the one-to-one correspondence between and )(rVext )(r guarantees that all ground state properties are determined by the electron density )(r In other words, )(r is the basic variable of the system at ground state and all physical quantities can be expressed as functionals of )(r The constant potential is trivial for ground state, since if 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 non-equilibrium 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 0V 0V )(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 non-equilibrium open system. Then we propose a modification of the Hohenberg-Kohn 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 non-equilibrium open systems.

PAGE 76

65 5.2 Hohenberg-Kohn 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. 1-1. 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. 1-1) 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 first-principles 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 L 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 )(r can uniquely determine the transport properties of the system, specifically, whether )(r alone is sufficient in determining the current through the system. Because the current is a function of the bias voltage, if )(r is the basic variable of the system in the same way it is for a ground state, knowing )(r should allow one to determine the voltage bias across the scattering region.

PAGE 77

66 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 Hohenberg-Kohn theorem for the ground state can apply at least in the two leads. Then, according to the theorem, the Fermi energy L is decided by the electron density in the left lead )(rL to an constant and the Fermi energy LV R is decided by the electron density in the right lead )(rR to another constant Since both leads are connected to a battery, these two constants, and should not be independent. However, according to the ground state Hohenberg-Kohn 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, and are independent if only given electron density. RV LV RV LV RV The result is that the voltage bias RL 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 )(r alone is not sufficient to describe all properties of the system, thus conventional DFT is not applicable. The Hohenberg-Kohn theorem must be modified in order to be applied to such systems. 5.2.1 Voltage-dependence of the Density Functional In order to see how the Hohenberg-Kohn theorem should be changed, we start from a rigorous approach to this problem based on the adiabatic time-dependent procedure. We take the same time-dependent process as that in reference 24 and reference 25: at time the unperturbed system consists of three uncoupled regions, the left lead, the t

PAGE 78

67 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 L and R 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 RL The couplings between leads and the scattering region turn on slowly from t to 0 t such that at any time, the system is in the stationary state. The goal is to relate the current at with the couplings completely turned on to the unperturbed Fermi energies 0t L and R Given the initial state, the properties of the time-dependent system can be obtained by solving the time-dependent Schrdinger equation )()()(ttHtti (5-1) The Hamiltonian can be written as ictijiijeeiiextiierVtVrVrVHtVHtH)1)((')()()(21)()(200 (5-2) where, in the Hamiltonian, the first term 0H consists of the kinetic energy, the nuclear-electron interaction and the electron-electron interaction The second term is an extra time-dependent external potential, where 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 the transmission probability of electrons tunneling through these two barriers is negligible; thus, the two leads and the middle barrier are well separated. The extV eeV 'V t

PAGE 79

68 constant c should be positive and small enough such that the couplings between leads and the middle barrier turn on adiabatically. For such a time-dependent system, Runge and Gross [69] showed that the wave function at time, t )(t can be determined by the time-dependent electron density )(t and the initial state, up to a time-dependent phase factor, )(),()()( tetti (5-3) where, (, is a functional of the electron density and the initial state. Particularly, if the initial state is the stationary ground state at t then according to the first Hohenberg-Kohn theorem for ground state, it is completely determined by the initial electron density; hence )( can be eliminated in the equation (5-3), so that )(t is determined up to a phase factor by only. Neglecting the phase factor, we come to the time-dependent DFT in this case, which claims that any time-dependent physical quantity can be expressed as the functional of time-dependent electron density. However, for our case, at t the system is not in the ground state. The initial state )( 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 R by the same amount, we know that the initial state is determined by the Hamiltonian and the voltage bias RLbV

PAGE 80

69 Applying the first Hohenberg-Kohn 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 and the voltage bias. bV At, the Equation (5-3) can be written as, 0t btiVtet),()()( (5-4) which means that at a given bV )(t is determined by to a time-dependent phase factor. bV, can be equivalently written as bVb| which means that the density functional form depends on the bias voltage We come to a conclusion that all steady-state properties at can be expressed as functionals of bV 0t at a given bias voltage. The most important of these functionals is the exchange-correlation energy functional bV xcE The voltage-dependence of those density functionals actually comes from the voltage-dependence of the initial state. The first Hohenberg-Kohn theorem for the equilibrium state is just a special case for 0 bV The voltage-dependence 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 voltage-dependence of the initial state varies with the system, the voltage-dependence 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

PAGE 81

70 the basic variable of the functional, the function )(r can not determine the bias voltage 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 bV can completely determine the system. 5.2.2 The Non-equilibrium Electron Density The key to eliminate the voltage-dependence of the functional in the equation (5-4) is to introduce other functions into the functional which can determine the bias voltage and thus determine the initial state bV )( together with the electron density For the non-equilibrium system, the electron density can be expressed as [20, 22], ),(21)(rGdir (5-5) where, is the non-equilibrium lesser Green's function defined for the stationary state. For the system under study, which is characterized by two Fermi energies ),(rG L and R (suppose RL ), the electron density can be divided into the equilibrium part eq and the non-equilibrium part n )()()(rrrneq (5-6) where, RRrGdirrGdirneq),(21)(),(21)( (5-7) The term eq is called the equilibrium density in the sense that it can be equivalently written as an equilibrium integral [20],

PAGE 82

71 RrGdrReq),(Im1)( (5-8) where, is the usual retarded Green's function. RG We want to show that if given both eq and n at the initial time, the bias voltage is determined and then the initial state is determined. When t the left lead is well separated with respect to other regions. The left lead is in a local equilibrium state with the Fermi energy L while the scattering region and the right lead are in a local equilibrium state with the Fermi energy R In this case, n in the right lead is zero and in the left lead it also can be written as an equilibrium integral, LRrGdrRn),(Im1)( (5-9) Since the left lead is in local equilibrium, the first Hohenberg-Kohn theorem for equilibrium state applies. We know that in the left lead, the Hamiltonian is determined by )()()(rrrneq to an arbitrary constant. If a constant is added to the Hamiltonian, changes to. Then, with 0V )(RG )(0VGR eq given, from equation (5-8), we see that R changes to 0VR Introduce the new R and into the equation (5-9), we see that with the given )(RG n L becomes 0VL Hence the bias voltage RLbV is completely determined without any arbitrary constant. Then, the Equation (5-4) can be rewritten as, )0(),0()0()0(ttetneqti (5-10)

PAGE 83

72 where, the voltage-dependence of the functional is eliminated. Based on the functional neq , all steady state properties can be expressed as functionals of eq and n 5.2.3 The Stationary Principle and the Single Particle Picture For the time-dependent system, we can define the action A as, 0)()()(dtttHtitA (5-11) where )(t is the many body wave function. If )(t is normalizable, then given the initial state, )( the stationary requirement of the action, 0)(tA (5-12) is completely equivalent to equation (5-3), where )(t is an arbitrary change of the wave function under the conditions that 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 (1-1). Thus, the Equation (5-12) can not be directly applied. However, for the case of time-dependent DFT, which means that the initial state is a stationary ground state, the action is the functional of the electron density AA As long as the wave function )(t yields the correct electron density, we still can use the Equation (5-12) to calculate the action and the stationary requirement becomes, 0)(tA (5-13) with the conditions that 0)0( and doesn't change the neutrality condition of two leads and the scattering region. This procedure can be generalized to the non-equilibrium

PAGE 84

73 initial state case. According to the previous subsection, for the open system under a finite bias, the action can be written as neqA and the stationary principle becomes, 0)()(tAtAneq (5-14) where the arbitrary change of electron densities should not change the initial conditions nor the neutrality condition of the system. For the adiabatic time-dependent process, the action reduces to 0)(dttEA (5-15) Then, the Equation (5-14) becomes, 0)()(tEtEneq (5-16) where, E is the total energy of the system at time t. Both the ground state and the time-dependent DFT work in the framework of a single-particle picture. We follow the same procedure. Consider an independent particle system, which has a same geometry as shown in Fig. 1-1 and whose single particle orbitals can be written as ,m where denotes all good quantum numbers other than the energy m At the initial time, since the two leads are well-separated with respect to the scattering region, ,m consists of three kinds of states, ,l ,r 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 time-dependent procedure described in previous sections will not change the occupation number of each state when it evolves with the time. Further

PAGE 85

74 assume that at the initial time, ,l ,r and ,b are perfectly normalized such that the neutrality conditions are satisfied in all three regions. Since the system is non-interacting, both Green's functions and are diagonalized by )(G )(rG ,m It is straightforward that equation (5-7) becomes, RLRmmmnmmeqf,2,,2,)( (5-17) where is the occupation number. Since the occupation number of each state does not change with time and mf RL n is, RLlln,2, (5-18) Assume that eq and n yielded by equaton (5-17) and (5-18) are the same as those of the true interacting system under study. Then, similar to the equilibrium state DFT, we write the total energy as VETE where T is the single-particle kinetic energy, RLRlllnmmmeqneqTTTTT,2,22121 (5-19) The potential energy can be written as, VE )()(121)()(212121rrrrrdrdEEErrVrdECoulxcCoulextV (5-20)

PAGE 86

75 where is the classical electron-electron repulsion energy and is the exchange-correlation energy. is the functional of CoulE xcE VE eq and n We can define two effective potentials based on as, VE nVneffeqVeqeffEVEV (5-21) The Equation (5-16) can be written in terms of those orbitals as, 0)21(0)21(,2,2RLRmneffmeqeffVV (5-22) There are two boundary conditions for the system under study. The open boundary condition incorporated with equation (5-22) 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 non-equilibrium electrons see a different mean field potential than the equilibrium electrons do. The self-consistent procedure makes two electron densities eq and n to converge respectively, and the current is carried by the non-equilibrium electrons. 5.3 Exchange energy for uniform electron gas The exchange-correlation energy can be written as the sum of the exchange part and the correlation part, cxxcEEE Both and are functionals of xE cE eq and n 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.

PAGE 87

76 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 )()(xlx then the single-particle orbitals can be written as, rkizkykxkizyxeVe l kkkzyx 2/1)(2/311),,( (5-23) where zyxzyxnlk,,,,2 with are integers. Set the axis to go from left to right, which is the direction of the current. There are two kinds of states for this system: the right-going states () and the left-going states ( zyxn,, z 0zk 0 zk ). From equations (5-17) and (5-18), we see that the occupation of those states is: below, all states are occupied and between and, only right-going states are occupied. Define two quantities and by ek ek tk ek tk LtRemkmk2,22222 Two electron densities eq (which comes from all states below ) and ek n (which comes from all right-going states above ), can be related to and After some algebra, we have ek ek tk 233 eeqk and )(61332etnkk We can see more clearly from this case that the total electron density neq alone is not enough to describe the system: two systems with the same have different electronic structures if they have different set of eq and n The exchange energy can be evaluated by the equation [30] xE 2122112),(141rdrdrrrEx (5-24) where ),(21rr is the density matrix and its diagonal part is the electron density. It can be

PAGE 88

77 calculated as krkieVrr 122),(21 where the summation is over all occupied For a large number of occupied states, the summation can be replaced by an integral, k 1212202/023200023214141),(rkikkrkikeddSindkkeddSindkkrrtee (5-25) where the first integral is the contribution of all states below which is an equilibrium state integral and can be evaluated analytically. We denote this term as ek eq The second integral is the contribution of occupied states between and, which can be denoted as ek tk n For the non-equilibrium integral, since only states with positive are occupied, the variable zk of the second integral goes from 0 to 2/ We introduce two variables 2/)(21rrR and Then the density matrix is only a function of and the exchange energy can be written as, 12rs s 2)(141ssKKRdEx (5-26) where K is a constant since the system is uniform. For a given electron density when the bias voltage increases, n increases and eq decreases until at a certain bias voltage maxV eq becomes zero and n equals to defines the biggest difference between maxV L and R we can have for a uniform system with a given at zero temperature. Given the electron density and the bias voltage RL we can numerically evaluate and K Fig. 5-1 shows the ratio as a function of total electron density and bias voltage. denotes the term calculated at zero bias. From the 0/KK 0K

PAGE 89

78 figure, we clearly see the voltage-dependence of the exchange energy for a given electron density This supports the analysis of the previous section: the electron density alone is not enough to describe the system. Also, the -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. =7.30 =33.77=17.29=2.16 Figure 5-1. Exchange energy at different bias voltages for different (Electrons/m 3 ) We notice that although the voltage-dependence 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. (5-26) as,

PAGE 90

79 sdsKsdsKsdsKKKKKneqneqneqnneqeqneqneq)(141141141**22 (5-27) The first term is the contribution of the exchange of equilibrium electrons. The second term comes from the non-equilibrium electrons and the third one denotes the contribution of exchanging between equilibrium and non-equilibrium electrons. The voltage-dependence of these three terms for eqK nK 29.7 Electrons/m 3 is shown in Fig. 5-2. At bias zero, since n is zero, there are no non-equilibrium electrons. Hence and are zero. As the bias voltage increases, nK neqK eq decreases, which leads to the decrease of, and eqK n increases, which leads to the increase of both and. The term reaches its maximum at a certain bias voltage and then decreases. At the maximum bias voltage, nK neqK neqK n and there are no equilibrium electrons, then both and are zero at this bias voltage. eqK neqK According to the previous section, the exchange energy can be expressed as a functional of and n We write K as, nfKK ,10 (5-28) where, 0K is the one for the equilibrium system with the electron density ; this term is well known. We plot the numerical results of as the function of 0/KK /n for different electron densities in Fig. 5-3.

PAGE 91

80 =7.30 (a) =7.30 (b) Figure 5-2. Split of the exchange energy: (a) the exchange between non-equilibrium electrons and the exchange between non-equilibrium and equilibrium electrons and (b) the exchange between equilibrium electrons.

PAGE 92

81 The numerical results show that we have a universal functional nf for all electron densities. These data can be fit to the function, nnbaf2 (5-29) where and 0497.0a 0436.0b =7.30 =2.16 =17.29 =33.77 / Figure 5-3. Exchange energy as the functional of /n for different electron densities. Then the exchange energy can be written as, nnxxbardCE223/41 (5-30) where, 7386.0xC Equation (5-30) gives the exchange energy functional of the uniform system for all bias voltages. It clearly shows the dependence on the non-equilibrium electron density as well as the dependence on the total electron density, which support the analysis of the

PAGE 93

82 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 eq and n 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 non-equilibrium 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.

PAGE 94

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 I-V curve can be compared to the experiments and peaks in the curve can be explained by tunneling through molecular orbitals. We apply the state-of-art approach which combines the non-equilibrium Greens functions technique and DFT to Au/S/Azobenzene/S/Au molecular junction. The results show that the azobenzene molecule is a good candidate of light-driven molecular switch, which is important for application of nanoelectronics. The I-V 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 bias-enhanced 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 non-equilibrium effects introduced by the 83

PAGE 95

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 voltage-dependence of exchange energy. By introducing non-equilibrium electron density, we show that the voltage-dependence 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 DE-FG02-02ER45995.

PAGE 96

APPENDIX A EXPANSION OF LESSER GREENS FUNCTION The non-equilibrium Greens functions are defined as )'1()1()'1,1(CTiG (A-1) where, the shorthand notation ),(111tx is used. The time-order operator is defined on two time branches, + branch, which runs from CT t to t and branch, which runs in the opposite direction. With two time labels, and which can be located on either of two time branches, we have four kinds of non-equilibrium Greens functions: 1t '1t )'1,1()'1,1()'1,1()'1,1()'1,1(~ttGGGGG '11'11'11'11,,,,tttttttt (A-2) These four Greens functions can be denoted as the time-ordered Greens function the greater Greens function the lesser Greens function and the antitime-ordered Greens function tG G G tG ~ For the Hamiltonian which is given by equation (2-1), a particular lesser Greens function defined as nkG, )()(,tcditGknnk (A-3) 85

PAGE 97

86 where, the indices ,k denote channels of the left lead and denotes the basis states of the device region. The angle brackets indicate the ensemble average. The time of the operator lies on the branch and the time of the operator lies on the + branch. n nd kC Figure A-1. Expansion of lesser Greens function The time of the perturbation can be either on + branch or on branch. nkG, The coupling term CouH is taken to be the perturbation V Then, the Dyson equation can be used to expand the Greens functions, where the term is the Greens function without the perturbation. Inserting the basis states into the Dyson equation, we have, VGGGG00 0G mmnmkkknkGVgnVGGknGknGkG,,00, (A-4) where, the term nGk0 is zero because electrons can not propagate from the left lead to the device region without the coupling term and the term is the bare Greens function of the left lead. The equation (A-4) can be used for both equilibrium and non-equilibrium Greens functions. For non-equilibrium case, since the time of the perturbation can be either on + branch or branch, the equation (A-4) gives two terms. kkg, + nt ktnt mmnmktkkGVg,,mtmnmkkkGVg ~ ,,ktVtVt +

PAGE 98

87 As an example, the expansion of the lesser Greens function contains two terms as shown in Fig. A-1. Thus, can be written as nkG, nkG, mtmnkkmntkkmknkGgGgVG][ ~ ,,,, (A-5)

PAGE 99

APPENDIX B DETAILS OF NON-EQUILIBRIUM SELF-CONSISTENT 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 self-energies of two leads according to equation (3-2). Since the tight-binding model is used to describe two leads, only the part L (R) of the device region is affected by the left (right) semi-infinite lead. Thus the dimension of the self-energy of the left (right) lead L ( 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, self-energies of two leads are not affected by the device region, so they dont change in the self-consistent procedure. The self-consistent procedure starts from an initial electron density )(0r in the device region. Given the electron density, the Hamiltonian can be calculated as )(][')'(')(21002rVVrrrrdrVHbiasxcion (B-1) where, the last term )(rVbias is the linear voltage drop, ydrVLRbias)( The current is flowing in y direction and d measures the distance between two leads. The retarded and advanced Greens function can be achieved by RLarHSG1, (B-2) where, S is the overlap matrix of Gaussian basis states and H is the Hamiltonian matrix. 88

PAGE 100

89 In order to see it more clearly, equation (B-2) can be expanded as 1RRRRRMRMRLRLRRMRMMMMMLMLMRLRLMLMLLLLLLrHSHSHSHSHSHSHSHSHSG (B-3) where, indices L, M, and R denote three parts of the device region respectively. Then, the lesser Greens function can be computed according to equation (2-11) and the density matrix is determined by dGiijij)(21 (B-4) where, indices j i, denote Gaussian basis states. The electron density can be directly related to the density matrix as jijijijirrrjjiirrrr,*,)()()( (B-5) where, )(ri is the Gaussian basis in real space. The electric static potential is defined as )(')'(')()(rVrrrrdrVrVbiasiones (B-6) An interesting quantity is the difference between the electric static potential under a finite bias and the electric static potential of zero bias, ,0esV ')'(')(0rrrrdrVVVVbiaseseses (B-7) which shows the effects of the voltage bias on the electric static potential.

PAGE 101

LIST OF REFERENCES [1] C. Joachim, J. Gimzewski, R. Schlittler, C. Chavy, Phy. Rev. Lett 74 2102 (1995). [2] M. A. Reed, C. Zhou, C. J. Muller, T. P. Burgin, and J. M. Tour, Science 278 252 (1997). [3] D. J. Wold, R. Haag, M. A. Rampi, and C. D. Frisbie, J. Phys. Chem. B 106 2813 (2002). [4] M. A. Reed, Proc. IEEE 97, 652 (1999). [5] C. Dekker, Phys. Today 22 (1999). [6] A. Aviram, and M. A. Ratner, Chem. Phys. Lett. 29 277 (1974). [7] N. D. Lang, Phys. Rev. B 52, 5335 (1995). [8] N. D. Lang, Phys. Rev. B 55, 4113 (1997). [9] N. D. Lang and Ph. Avouris, Phys. Rev. Lett. 81, 3515 (1998). [10] M. Di Ventra, S. T. Pantelides, and N. D. Lang, Phys. Rev. Lett. 84, 979 (2000). [11] N. D. Lang and Ph. Avouris, Phys. Rev. Lett. 84, 358 (2000). [12] C. C. Wan, J.-L. Mozos, G. Taraschi, J. Wang, and H. Guo, Appl. Phys. Lett. 71, 419 (1997). [13] J. Wang, H. Guo, J.-L. Mozos, C. C. Wan, G. Tarachi, and Q. Zheng, Phys. Rev. Lett. 80, 4277 (1998). [14] G. Tarachi, J.-L. Mozos, C. C. Wan, H. Guo, and J. Wang, Phys. Rev. B 58, 13138 (1998). [15] Yongqiang Xue, Supriyo Datta and Mark A. Ratner, Chem. Phys. 281, 151 (2000). [16] H. J. Choi, J. Ihm, Y. Y. Yoon, and S. G. Louie, Phys. Rev. B 60, 14009 (1999). 90

PAGE 102

91 [17] U. Landman, R. N. Barnett, A. G. Scherbakoy, and P. Avouris, Phys. Rev. Lett. 85, 1958 (2000). [18] B. Larade, J. Taylor, H. Mehrez, and H. Guo, Phys. Rev. B 64, 075420 (2001). [19] H. J. Choi, J. Ihm, S. G. Louie, and M. L. Cohen, Phys. Rev. Lett. 84, 2917 (2000). [20] J. Taylor, H. Guo, and J. Wang, Phys. Rev. B 63, 121104 (2001). [21] H. Mehrez, Alex. Wlasenko, Brian Larade, Jeremy Taylor, Peter Grutter, and Hong Guo, Phys. Rev. B 65, 195419 (2002). [22] Y. Xue, S. Datta, S. Hong, R. Reifenberger, J. I. Henderson, and C. P. Kubiak, Phys. Rev. B 59, 7852 (1999). [23] P. Hohenberg, W. Kohn, Phys. Rev. 136, B864 (1964). [24] Yigal Meir and Ned S. Wingreen, Phys. Rev. Lett. 68, 2512 (1992). [25] C. Caroli, R. Combescot, Pl. Nozieres, and D. Saint-James, J. Phys. C 4, 916 (1971). [26] C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1998). [27] A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Publications, New York, 1996). [28] D. E. Bernholdt, E. Apra, H. A. Fruchtl, M. F. Guest, R. J. Harrison, R. A. Kendall, R. A. Kutteh, X. Long, J. B. Nicholas, J. A. Nichols, H. L. Taylor, A. T. Wong, G. I. Fann, R. J. Littlefield, and J. Nieplocha, Int. J. Quantum Chem., Quantum Chem. Symp. 29, 475-483 (1995). [29] R. M. Dreizler and E. K. U. Gross, Density Functional Theory (Springer-Verlag, Berlin, Germany, 1990). [30] R. G. Parr and W. T. Yang, Density-Functional Theory of Atoms and Molecules (Oxford University Press, Oxford, UK, 1989). [31] Francesca Moresco, Gerhard Meyer, and Karl-Heinz Rieder, Phys. Rev. Lett.. 86, 672 (2000). [32] J. M. Seminario, P. A. Derosa, and J. L. Bastos, J. Am. Chem. Soc. 124, 10266 (2001). [33] R. Gutierrez, G. Fagas, G. Cuniberti, F. Grossmann, R. Schmidt, K. Richter, Phys. Rev. B 65, 113410 (2002). [34] P. Orellana and F. Claro, Phys. Rev. Lett. 90, 178302 (2003).

PAGE 103

92 [35] E. G. Emberly and G. Kirczenow, Physl. Rev. Lett. 91, 188301 (2003). [36] S. Monti, G. Orlandi, and P. Palmieri, Chem. Phys. 71,87 (1982). [37] H. Rau, J. Photochem. 26, 221 (1984). [38] T. Nagele, R. Hoche, W. Zinth, and H. Wachtveitl, Chem. Phys. Lett. 272, 489 (1997). [39] F. W. Schulze, H. J. Petrick, H. K. Cammenga, and H. Klinge, Z. Phys. Chem. Neue Fol. 107, 1 (1977). [40] E. R. Talaty and J. C. Fargo, Chem. Commun. 2, 65 (1967). [41] T. Hugel, N. Holland, A. Cattani, L. Moroder, M. Seitz, H. Gaub, Science 296, 1103 (2002). [42] P. S. Krstic, X.-G. Zhang, and W. H. Buttler, Phys. Rev. B 66, 205319 (2002). [43] S. Sanvito, C. J. Lambert, J. H. Jefferson, A. M. Bratkovsky, Phys. Rev. B 59, 11936 (1999). [44] R. Landauer, IBM J. Res. Dev. 1, 233 (1957). [45] S. A. Rishton, Y. Lu, R. A. Altman, A. C. Marley, X. P. Bian, C. Jahnes, R. Viswannathan, G. Xiao, W. J. Gallagher, and S. S. P. Parkin, Microelectron. Eng. 35, 249 (1997). [46] S. S. Parkin, K. P. Roche, M. G. Samant, P. M. Rise, R. B. Beyers, R. E. Scheuerlein, E. J. O'Sullivan, S. L. Brown, J. Bucchigano, D. W. Abraham, Y. Lu, M. Rooks, Pl L. Trouilloud, R. A. Wanner, and W. J. Gallagher, J. Appl. Phys. 85, 5828 (1999). [47] J. S. Moodera, E. F. Gallagher, K. Robinson, and J. Nowak, Appl. Phys. Lett. 70, 3050 (1997). [48] C. T. Tanaka, J. Nowak, and J. S. Moodera, J. Appl. Phys. 86, 6239 (1999). [49] J. Zhang and R. M. White, J. Appl. Phys. 83, 6512 (1998). [50] P. Zahn and I. Mertig, Phys. Rev. Lett. 75, 2996 (1995). [51] J. M. MacLaren, X.-G. Zhang, W. H. Butler, and Xindong Wang, Phys. Rev. B 59, 5470 (1999). [52] W. H. Butler, X.-G. Zhang, T. C. Schulthess, and J. M. MacLaren, Phys. Rev. B 63, 054416 (2001). [53] J. Mathon and A. Umerski, Phys. Rev. B 63, 220403(R) (2001).

PAGE 104

93 [54] O. Wunnicke, N. Papanikolaou, R. Zeller, P. H. Dederichs, V. Drchal, and J. Kudmovski, Phys. Rev. B 65, 064425 (2002). [55] Xindong Wang, J. Appl. Phys. 83, 6518 (1998). [56] A. H. Davis and J. M. MacLaren, J. Appl. Phys. 87, 5224 (2000). [57] T. Nagahama, S. Yuasa, Y. Suzuki, and E. Tamura, J. Appl. Phys. 91, 7035 (2002). [58] S. Zhang, P. M. Levy, A. C. Marley, and S. S. P. Parkin, Phys. Rev. Lett. 79, 3744 (1997). [59] J. M. MacLaren, X.-G. Zhang, and W. H. Butler, Phys. Rev. B 56, 11827 (1997). [60] W. H. Butler, X.-G. Zhang, T. C. Schulthess, and J. M. MacLaren, Phys. Rev. B 63, 092402 (2001). [61] J. M. MacLaren, S. Crampin, D. D. Vvednsky, R. C. Albers, and J. B. Pendry, Comput. Phys. Commun. 60, 365 (1990). [62] W. Wulfhekel, M. Klaua, D. Ullmann, F. Zavaliche, J. Kirschner, R. Urban, R. T. Monchesky, and B. Heinrich, Appl. Phys. Lett. 78, 509 (2001). [63] H. L. Meyerheim, R. Popescu, J. Kirschner, N. Jedrecy, M. Sauvage-Simkin, B. Heinrich, and R. Pinchaux, Phys. Rev. Lett. 87, 076102 (2001). [64] H. L. Meyerheim, R. Popescu, N. Jedrecy, M. Sauvage-Simkin, R. Pinchaux, B. Heinrich, and J. Kirschner, Phys. Rev. B 65, 144433 (2002). [65] M. Bowen, V. Cros, F. Petroff, A. Fert, C. M. Boubeta, J. L. Costa-Kramer, J. V. Anguita, A. Cebollada, F. Briones, J. M. DeTeresa, L. Morellon, M. R. Ibarra, F. Guell, F. Peiro, and A. Cornet, Appl. Phys. Lett. 79, 1655 (2001). [66] X.-G. Zhang, W. H. Butler, and A. Bandyopadhyay, Phys. Rev. B 68, 092402 (2003). [67] M. Brandbyge, J. Mozos, P. Ordejon, J. Taylor, K. Stokbro, Phys. Rev. B 65, 165401 (2002). [68] T. Ando, Y.Arakawa, K. Furuya, S. Komiyama and H. Nakashima, Mesoscopic Physics and Electronics (Springer-Verlag, Berlin, 1998). [69] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).

PAGE 105

BIOGRAPHICAL SKETCH Chun Zhang was born on August 3, 1973, in Yongshun County, Hunan Province, China. He got his B. S degree in Nuclear Science in Fudan University, Shanghai, China in 1996 and M. S degree in physics in the same University in 2000. He entered the physics department of University of Florida in August, 2000. 94