Table Of Content

Conductance calculations for quantum wires and interfaces: mode matching and Green functions P. A. Khomyakov, G. Brocks,∗ V. Karpan, M. Zwierzycki, and P. J. Kelly Computational Materials Science, Faculty of Science and Technology and MESA+ Research Institute, University of Twente, P. O. Box 217, 7500 AE Enschede, The Netherlands. 5 (Dated: February 2, 2008) 0 Landauer’s formularelates theconductanceofaquantumwire orinterfacetotransmission prob- 0 abilities. Totaltransmission probabilitiesarefrequentlycalculated usingGreenfunctiontechniques 2 and an expression first derived by Caroli. Alternatively, partial transmission probabilities can be n calculated from the scattering wave functions that are obtained by matching thewave functions in a the scattering region to the Bloch modes of ideal bulk leads. An elegant technique for doing this, J formulated originally by Ando, is here generalized to any Hamiltonian that can be represented in 5 tight-bindingform. Amorecompactexpressionforthetransmissionmatrixelementsisderivedand 2 it is shown how all the Green function results can be derived from the mode matching technique. We illustrate this for a simple model which can be studied analytically, and for an Fe|vacuum|Fe ] tunneljunction which we study using first-principlescalculations. i c s PACSnumbers: 73.63.-b,73.40.-c,73.20.-r,85.35.-p - l r t m I. INTRODUCTION matching17 and Green function15,16,22 approaches is not immediately obvious. Indeed, it wasrecently statedthat . at Since the discovery of giant magnetoresistance in Ando’s approach is incomplete and does not yield the m metallic multilayers there has been considerable interest correct expression for the conductance.22 - in studying electronic transport in layered materials.1,2 d At the same time, experimental control of the lateral n scalehasenabledstudies ofelectronictransportinquan- Inthis paper we demonstrate that the twoapproaches o tum wires of atomic dimensions.3 Because of the small are completely equivalent. In the Green function ap- c [ dimensions involved, the transport properties of both of proach,asmallimaginarypartmustbe addedtoorsub- these systems should be understoodon the basis oftheir tracted from the energy in order to distinguish between 1 atomic structure. This perception has generated a large the retarded and advanced forms.5,6,7,10,11,12,13,14,22 In v effortinrecentyearstocalculatetheconductanceofmul- mode matching,scatteringwavefunctions arecalculated 9 0 tilayers and quantum wires from first principles. Several thatincorporatetheretardedoradvancedboundarycon- 6 different approaches have been formulated which have a ditions directly. This makes it possible to solve the 1 common basis in the Landauer-Bu¨ttiker approachor are scattering problem also at real, instead of complex en- 0 equivalent to it. In the linear response regime, the con- ergies. In addition to yielding the total conductance, by 5 ductance G is expressed as a quantum mechanical scat- focussingonwavefunctionsthecontributionofeachindi- 0 tering problem4 and can be simply related to the total vidualscatteringchannelcanbeidentified. Inparticular, / t transmissionprobabilityatthe Fermienergy,T(E ), as we derive a simple, compactexpressionfor the transmis- a F m sion matrix elements, see Eq. (67). e2 G = T(E ). (1) - F h d n The multilayer or quantum wire is generally considered The paper is organized as follows. In the next section o as a scattering region of finite size, sandwiched between the Hamiltonian we will use is introduced. This model c two semi-infinite ballistic wires. Aiming at a materials allows us to study both quantum wires that are finite : v specific description, most current approaches treat the in the directions perpendicular to the wire, and systems i X electronicstructurewithintheframeworkofdensityfunc- that are periodic in these directions such as single inter- tional theory (DFT).5,6,7,8,9,10,11,12,13,14 faces,sandwichesandmultilayers. We willuse the single r a FrequentlytheconductanceiscalculatedusingaGreen term “quantum wire” to describe both systems. In Secs. function expression first derived by Caroli et al.15,16 III and IV the mode matching and Green function tech- An alternative technique, suitable for Hamiltonians that niques are summarized. The equivalence of the trans- can be represented in tight-binding form, has been for- mission matrices obtained using these two techniques is mulated by Ando.17 It is based upon directly match- demonstratedinSec. VandtheCaroliexpressionforthe ing the wave function in the scattering region to the conductance is derived from the mode matching expres- Bloch modes of the leads. The latter technique has sions. InSec. VIthetwotechniquesareappliedfirsttoa been applied to conductance calculations at the empiri- simple analytical model,23 and then to an Fe|vacuum|Fe caltight-binding level,18 as wellas onthe first-principles tunnel junction using numerical first-principles calcula- DFT level.8,19,20,21 The relationship between the mode tions. ThemainconclusionsaresummarizedinSec. VII. 2 II. HAMILTONIAN Wesetupatight-bindingrepresentationofthe Hamil- tonian. This is not a severe restriction since a first- principles DFT implementation that uses a representa- tion on an atomic orbital basis set has the same math- ematical structure as a tight-binding model.24,25 Alter- natively, an implementation that uses a representation of the Hamiltonian on a grid in real space, can also be mappedontoatight-bindingmodel.26Webeginbydivid- FIG. 1: (Color online) Hamilton matrix of a quantum wire divided into slices. The left (L) and right (R) leads are ideal ing the system into slices (“principal layers”)perpendic- periodic wires that span the cells i = −∞,...,0 and i = ular to the wire direction.27 The thickness of these slices S+1,...,∞, respectively. The scattering region spans cells is chosen such that there is only an interaction between i=1,...,S. neighboring slices. Labelling each slice with an index i, the Schr¨odinger equation of the quantum wire becomes A. Bloch matrices −H c +(EI−H )c −H c =0, (2) i,i−1 i−1 i,i i i,i+1 i+1 Thefirststepconsistsoffindingsolutionsfortheleads, fori=−∞,...,∞. AssumingthateachslicecontainsN for which Eq. (2) can be simplified to sites and/or orbitals, c is a vector of dimension N con- i taining the wave function coefficients on all sites and/or −Bc +(EI−H)c −B†c =0. (3) i−1 i i+1 orbitals of slice i. The N ×N matrices H and H i,i i,i±1 consist of on-slice and hopping matrix elements of the These equations hold for i = −∞,...,−1 and i = Hamiltonian, respectively. I is the N ×N identity ma- S +2,...,∞, i.e. the left and right leads. To keep the trix. A schematic representation of the structure of the notation as simple as possible, we have omitted the sub- Hamiltonian is given in Fig. 1. scripts L and R, see Fig. 1. Since the leads are peri- Eq. (2) is valid both for quantum wires that are finite odic wires, one can make the ansatz that the solutions in the directions parallel to the slices, and for layered have Bloch symmetry, i.e. ci = λci−1, ci+1 = λ2ci−1, systems that are periodic in these directions. In the lat- whereλistheBlochfactor. SubstitutingthisintoEq.(3) ter case, translations in the transverse direction can be results in a quadratic eigenvalue equation of dimension described in terms of a Bloch wave vector in the two- N. The latter canbe solvedmost easilyby transforming dimensionalBrillouinzone,k ,whichisagoodquantum it to an equivalent linear (generalized) eigenvalue prob- k number and the system becomes effectively one dimen- lem of dimension 2N and solving this by a standard sional. Explicit expressions for the Hamilton matrix el- algorithm.29,30 ements depend upon the particular localized orbital ba- It can be shown that the equation generally has 2N sis or realspace grid representationused.28 Since details solutions,whichcanbedividedintoN right-goingmodes of the tight-binding muffin tin orbital scheme used in and N left-going modes,31 labelled as (+) and (−) in Refs. 8,19,20,21 are given in Ref. 25 and of the real- the following. Right-going modes are either evanescent space high-order finite difference method can be found wavesthataredecayingtotheright,orwavesofconstant in Ref. 26, they will not be discussed further here. amplitudethatarepropagatingtotheright,whereasleft- The system is divided into three parts, with i = goingmodes aredecaying or propagatingto the left. We −∞,...,0correspondingtotheleftlead(L),i=1,...,S denote the eigenvalues by λn(±) where n = 1,...,N, to the scattering region (S) and i = S + 1,...,∞ to the corresponding eigenvectors by un(±) and write the the right lead (R). The leads are assumed to be ideal eigenvalue equation as wires characterized by a periodic potential. It is then −Bu (±)+(EI−H)λ (±)u (±)−B†λ (±)2u (±)=0. n n n n n sufficient to identify a slice with a translational period (4) along the wire. By construction, the Hamilton matrix is In the following we assume that the vectors u (±) are the same for each slice of the leads, i.e. H ≡ H , n i,i L/R normalized;notethatingeneraltheyarenot orthogonal. Hi,i−1 ≡ BL/R and Hi,i+1 ≡ B†L/R for the left/right One can distinguish right- from left-going evanescent leads. Fig. 1 summarizes our model of a quantum wire. modes on the basis of their eigenvalues; right going evanescentmodes have |λ(+)|<1 and left going evanes- cent modes have |λ(−)| > 1. Propagating modes have III. MODE MATCHING APPROACH |λ(±)| = 1, so here one has to determine the Bloch ve- locityanduseitssigntodistinguishrightfromleftprop- Eq. (2) can be solvedin a particularly convenient way agation. TheBlochvelocitiesaregivenbytheexpression by a technique which we call mode matching. In this, we follow Ando and generalize his approach to a slice 2a geometry.17 vn(±)=− ~ Im λn(±)un(±)†B†un(±) , (5) (cid:2) (cid:3) 3 where a is the slice thickness. A derivation of this equa- with a modified Hamilton matrix defined so H′ =H i,j i,j tion is given in Appendix A. for the elements {i,j =0,1}, {i=1,...,S; j =i,i±1} Since the eigenvectors are non-orthogonal,it is conve- and {i,j =S+1,S}, but with nient to define dual vectors u (±) n H′ = H +B F−1(−), u†(±)u (±)=δ ; u†(±)u (±)=δ . (6) 0,0 L L L n m n,m e n m n,m H′ = H +B†F (+). (13) S+1,S+1 R R R Any wave function in the leads can be expressed as a e e linear combination of the lead modes. This can be done H′ =0 andH′ =0,so the modifiedscattering 0,−1 S+1,S+2 in a compact form using the two N ×N Bloch matrices regionisdecoupledfromtheleads. Ontheleft-handside for right- and left-going solutions of Eq. (12), we have a source term with N Q =B F−1(+)−F−1(−) , (14) F(±)= λ (±)u (±)u†(±). (7) 0 L L L n n n n=1 (cid:2) (cid:3) X and Q = 0 for i = 1,...,S +1. Eq. (12) defines a set i For any integer i, Fi is givenby a simeilar expressionbut of (S +2)×N linear equations. Because the Hamilton with λ replaced by λi. From the foregoing it is easy to matrixisblocktridiagonal,eachblockbeingofdimension n n show that the Bloch matrices obey the equation N, this set of equations can be solved efficiently using a block Gaussian elimination scheme.30 The total wave −BF−1(±)+(EI−H)−B†F(±)=0. (8) function c can then be obtained by back substitution. i The transmission is obtained from the wave function A general solution in the leads can be expressed in in the right lead c (+). In particular, choosing the terms of a recursion relation S+1 incoming wave as one of the modes of the left lead, i.e. ci = ci(+)+ ci(−) c0(+) = uL,m(+), generalized transmission matrix ele- = Fi−j(+)c (+)+Fi−j(−)c (−). (9) ments τn,m are defined by expanding cS+1(+) in modes j j of the right lead Fixingthecoefficientsinoneslicethensetstheboundary conditionsanddeterminesthesolutionforthewholelead. N c (+)= u (+)τ . (15) S+1 R,n n,m n=1 X B. Transmission matrix Byletting c (+) runoverallpossible incomingmodes of 0 the left lead u (+); m = 1,...,N, a full transmission The scattering region is defined by i = 1,...,S, see L,m matrix is obtained. Fig. 1. Right and left of the scattering regions one has Matrix elements can be defined for all modes, prop- the recursion relations for the states in the leads from agating and evanescent, but of course only matrix ele- Eq. (9) ments where n,m denote propagating modes contribute c = F−1(+)c (+)+F−1(−)c (−) to the real physical transmission. These modes can be −1 L 0 L 0 = F−1(+)−F−1(−) c (+)+F−1(−)c ,(10) selected by making use of their eigenvalues; see the dis- L L 0 L 0 cussionfollowingEq.(4). Thephysicaltransmissionma- with c = c(cid:2)(+)+ c (−) and (cid:3) trixelementsarethenfoundbynormalizingwithrespect 0 0 0 to the current16 c =F (+)c (+)+F (−)c (−), (11) S+2 R S+1 R S+1 v (+)a where the subscripts L and R distinguish between the t = R,n L τ , (16) n,m n,m Bloch matrices of the left and right leads. svL,m(+)aR The boundary conditions for the scattering problem are set up in the usual way. The vector c0(+) is treated where vL,m(+) and vR,n(+) are the Bloch velocities in as the source, i.e. as a fixed incoming wave from the left the directionofthe wireforthe right-propagatingmodes lead. There is no incoming wave from the right lead, so m and n in the left and right leads, respectively, see we set cS+1(−)=0. Eq. (5); aL and aR are the slice thicknesses of left and Havingsettheboundaryconditions,Eqs.(10)and(11) right leads.32 The total transmissionprobability is given can be used to rewrite Eq. (2) in the region not covered by by Eq. (3), i.e. for i = 0,...,S +1. This describes the wave function in the scattering region and the matching (+) to the solutions in the leads. Eq. (2) in this region is T(E)= |tn,m|2, (17) rewritten as n,m X −H′ c + EI−H′ c −H′ c =Q c (+), and the conductance is given by Eq. (1) evaluated at i,i−1 i−1 i,i i i,i+1 i+1 i 0 (12) E =E . F (cid:0) (cid:1) 4 C. Green function matrix Here g (z) and g (z) are the surface Green functions of L R the semi-infinite left and right leads, respectively, which can be calculated using an iterative technique: denoting Solving the set of linear equations Eq. (12) directly leads to the conductance. However, to facilitate a con- G[n](z)asthesolutionofanequationsimilartoEq.(21), i,j nection to the Green function approach discussed in but with H = 0 for {i > n∨j > n}, one can easily i,j the next section, we can formulate the solution in a derive the right-going recursion relation slightly different way. A finite Green function matrix G′ (z), i,j =0,...,S+1 can be defined by zI−Hn+1,n+1−Hn+1,nG[nn,]n(z)Hn,n+1 i,j h G[n+1] (z)=Ii. (23) −H′ G′ + zI−H′ G′ −H′ G′ =Iδ , n+1,n+1 i,i−1 i−1,j i,i i,j i,i+1 i+1,j i,j (18) Foranidealwirewithi,j =−∞,...,n,G[n] (z)=g (z) with z complex. (cid:0)Note, how(cid:1)ever, that the matrices H′0,0 should be independent of n resulting inn,tnhe followLing and H′S+1,S+1 are non-Hermitian and G′i,j(E) is also equation for the surface Green function, uniquely defined for real energies. The Green function matrix allows the solution of Eq. (12) to be written as zI−HL−BLgL(z)B†L gL(z)=I. (24) h i c =G′ (E)Q c (+). (19) Several iterative algorithms have been formulated for i i,0 0 0 solving this non-linear matrix equation.27,35,36 A simi- Asbefore,thetransmissioncanbeextractedati=S+1 lar reasoning based upon a left-going recursion for the and comparison with Eq. (15) gives right lead results in an equation for the surface Green function g (z) of the right lead R τn,m =u†R,n(+)G′S+1,0(E)Q0uL,m(+). (20) zI−H −B†g (z)B g (z)=I. (25) R R R R R which can be useed in Eq. (16). The Green function ma- Again, sethting z = E +iη in Eqs.i(24) and (25) defines trix block G′ (E) can be found by solving Eq. (18) the usualretardedsurface Greenfunctions g (E). Al- S+1,0 L/R using a recursive algorithm that resembles a Gaussian though we are mainly interested in the physical limit elimination scheme.25,33 lim , in practice a finite value of η is often used in η→0 order to make the iterative algorithms stable. The quantities IV. GREEN FUNCTION APPROACH Σ (E)=B g (E)B† ; Σ (E)=B†g (E)B , (26) L L L L R R R R An apparently quite different route to the transmis- which appear in Eqs. (22)-(25), are called the self- sion matrix starts by defining an infinite Green function energies of the left and right leads, respectively.16 Once matrix G (z) for i,j =−∞,...,∞ with respect to the theseareobtained,thefiniteHamiltonmatrixofEq.(22) i,j original Hamiltonian of Eq. (2). is constructed and the retarded Green function matrix G (E) can be found using a recursive algorithm.33 Us- i,j −H G +(zI−H )G −H G =Iδ . ingtheleadmodesthe transmissionmatrixelementscan i,i−1 i−1,j i,i i,j i,i+1 i+1,j i,j (21) then be calculated, as will be shown in the next section, Choosing z = lim (E ±iη) defines as usual the re- Sec.IVB. Alternatively,thetotaltransmissionprobabil- η→0 tarded/advanced Green function matrix. We shall use ity can be expressed in a form that does not require the G (E) to denote the retarded Green function matrix lead modes explicitly, which is discussed in Sec. V. i,j and Ga (E) to denote the advanced Green function ma- i,j trix. B. Transmission matrix A. Partitioning The transmission matrix can be obtained from the Green function matrix of Eq. (21). To do this, we adaptaFisher-Leetypeofapproachtoourtight-binding Eq.(21)ismostconvenientlysolvedbyapplyingapar- formulation.37 Assuming that the unperturbed reference titioningtechnique.16,34Itisstraightforwardtoshowthat wave function is the Bloch mode u (+) that comes in thefinitepartG (z)definedfori,j =0,...,S+1canbe L,m i,j from the left lead, the Lippmann-Schwinger equation38 derivedfromaclosedsetofequations,similartoEq.(21), in tight-binding form is butwithH replacedbyH′′ whereH′′ =H forthe i,j i,j i,j i,j elements {i,j = 0,1}, {i = 1,...,S; j = i,i±1} and c = u (+)+ G V u (+) i L,m,i i,j j,k L,m,k {i,j =S+1,S}, but with j,k X H′′ (z) = H +B g (z)B†, 0,0 L L L L = Fi(+)+ G V Fk(+) u (+).(27) H′S′+1,S+1(z) = HR+B†RgR(z)BR. (22)  L Xj,k i,j j,k L  L,m   5 Hereu (+)isthereferencewavefunctioninslicei. It terms of its eigenmodes. The columns of such a Green L,m,i obeys Bloch symmetry and u (+) ≡ u (+) is the function obey the equation L,m L,m,0 Blochmodeattheorigin,seeSec.IIIA. ThematrixV j,k represents the perturbationwith respect to the ideal left −BG(0) + E+I−H G(0)−B†G(0) =Iδ , (31) i−1,j i,j i+1,j i,j lead. Eq. (27) can be simplified using the Dyson equation, where E+ = E(cid:0) + iη. F(cid:1)or i 6= j the solution is simi- which in tight-binding form reads lar to that of the wave functions, see Eq. (3). In ad- dition, the retarded Green function should consist only G = G(0)+ G V G(0) i,0 i,0 i,j j,k k,0 of propagating waves that move outwards from the δ- Xj,k source and/or evanescent states that decay away from the source.38 From Eq. (9), we have the column recur- = Fi(+)+ G V Fk(+) G(0), (28) sion relations  L i,j j,k L  0,0 j,k X G(0)(E) = Fi−j(−)G(0)(E), i<j,   i,j j,j using Eq. (32). Comparing Eqs. (27) and (28) one finds the simple expression G(i,0j)(E) = Fi−j(+)G(j0,j)(E), i>j. (32) −1 ci =Gi,0(E) G(00,0)(E) uL,m(+). (29) ThediagonalblockG(j0,j)(E)cannowbeobtainedbycom- h i bining Eqs. (31) and (32), which gives for i=j Fromthedefinitionofthegeneralizedtransmissionma- trix elements, cf. Eq. (15), one then obtains the expres- −BF−1(−)+E+I−H−B†F(+) G(0) =I. (33) sion j,j τ =u† (+)G (E) G(0)(E) −1u (+). (30) Elim(cid:2)inating E+I−H using Eq. (8) th(cid:3)en yields n,m R,n S+1,0 0,0 L,m h i −1 To find τne,m one needs to calculate only the Green func- G(j0,j)(E) =B F−1(+)−F−1(−) , (34) tion matrix blocks G (E) of the full system and S+1,0 h i h i G(0)(E) of the idealleft lead. The physicaltransmission or the equivalent 0,0 matrix elements and the total transmission probability can then be obtained from Eqs. (16) and (17). G(0)(E) −1 =B†[F(−)−F(+)], (35) j,j h i Eqs. (32) and (34) represent the full expression for the V. MODE MATCHING VERSUS GREEN FUNCTIONS Green function G(0) of an infinite ideal wire in terms i,j of the Bloch matrices F(±) and thus in terms of the The two seemingly different formalisms introduced in eigenmodes. Note that we can set E+ = E since, in Secs. IIIandIVareinfactcloselyrelated. Inthissection terms of the modes, the retarded Green function matrix we will show how all Green function results can be ob- is uniquely defined for real energies. tained from the mode matching approach. We begin by The advanced Green function matrix G(0)a(E) can be i,0 expressing the Green function matrices of ideal wires in found from a similar procedure. It should consist of termsoftheBlochmatrices,F(±). Theseexpressionsare propagatingwaves that move towards the source and/or thenusedtoprovethatthetransmissionmatrixelements evanescentstatesthatgrowtowardsthesource. Onecan obtainedfromthemodematchingandGreenfunctionap- constructtwonewBlochmatricesFa(±),whicharesim- proaches,cf. Eqs.(20)and(30),areinfactidentical. Af- ilar to those defined in Eq. (7). In Fa(+) one collects ter that, we show that the transmission matrix elements the modes thatare decayingto the right(growingto the are independent of the exact positions within the leads left) and modes that are propagating to the left. Fa(−) thatareusedtomatchtheleadstothescatteringregion, thencontainsmodesthatgroworpropagatetotheright. apart from a trivial phase factor. Then we derive from the mode matching expressionfor the total transmission N Fa(±) = λa(±)ua(±)ua†(±), with (36) probability the Green function expression known as the n n n Caroli expression.15 Finally, a more compact expression nX=1 for the transmission matrix elements is derived. λan(±) = λn(∓), uan(±)=eun(∓)propagating λa(±) = λ (±), ua(±)=u (±)evanescent. n n n n A. Green functions of ideal wires in terms of Bloch Using these definitions, expressions for the advanced matrices Green function matrix are obtained from Eqs. (32) and (34) by replacing F(±) with Fa(±). We begin by deriving an expression for the retarded From the general relation between retarded and ad- Green function matrix G(0) of an ideal infinite wire in vanced Green functions, G = Ga †, the following i,j i,j j,i (cid:0) (cid:1) 6 row recursion relations can be deduced for the retarded Inotherwords,thetwoGreenfunctionsdiscussedinSecs. Green function III and IV are identical. By comparing Eqs. (14) and (34) one has G(0)(E) = G(0)(E) Fa†(+) j−i, i<j, i,j i,i −1 G(0)(E) = G(0)(E)(cid:2)Fa†(−)(cid:3)j−i, i>j. (37) Q0 = G(00,0)(E) . (46) i,j i,i h i TheretardedGreenfunctio(cid:2)nG(s)(E(cid:3))ofasemi-infinite In conclusion, the two expressions for the generalized i,0 transmission matrix elements, Eqs. (20) and (30), are wire extendingfromi=−∞,...,0canbeobtainedusing identical. a similar technique. Instead of Eq. (32), we get G(s)(E)=Fi(−)g(E), i<0. (38) i,0 C. Invariance of transmission probability where g(E) = G(s)(E) is the surface Green function. 0,0 Apartfromtrivialphase factors,the transmissionma- Using this in Eq. (31) gives for i = 0 and j = 0 and for trix elements should not depend on where in the ideal i=−1 and j =0, respectively, lead the wave function matching is carried out. In a re- centpaperit wasstatedthat Ando’sexpressionfort , −BF−1(−)+E+I−H g=I, n,m Eqs. (20) and (16), lacks this invariance property and is (cid:2)−BF−1(−)+E+I−H(cid:3)F−1(−)g =B†g. (39) therefore incomplete.22 One can however prove directly from Eq. (20) that the transmission matrix elements do Note tha(cid:2)t the B† term is absen(cid:3)t in the first equation have the required invariance property.39 The proof be- since we aredealing with a semi-infinite wire. These two comes easier if the equivalence of Eqs. (20) and (30), equations can be easily solved to find an expression for established above, is used. the surface Green function The scattering region runs from slices 0 to S + 1 if we include the boundaries with the left and right leads. g(E)=F−1(−) B† −1. (40) This means that the Green function matrix G with i,j indicesi,j outsidethisregionobeystheequationsforthe Eqs.(38)and(40)representtheG(cid:0)ree(cid:1)nfunctionofasemi- idealleads. Fromthecolumnandrowrecursionrelations, infinite ideal wire extending from i = −∞,...,0. In Eqs. (32) and (37), one derives a similar fashion, one gets for the Green function of a semi-infinite ideal wire extending from i=0,...,∞ G (E)=Fi (+)G (E) Fa†(−) j, (47) S+1+i,j R S+1,0 L G(s)(E)=Fi+1(+)B−1, i≥0. (41) h i i,0 for j < 0, i > 0. In a similar way, one derives for the Green function matrix of the left lead Analogously to Eqs. (38)-(41), one can define the ad- vanced Green function matrix Gi(,s0)a(E) in terms of the G(0)(E) = Fj(+)G(0)(E) Fa†(−) j, (48) j,j L 0,0 L Bloch matrices Fa(±). Moreover, since [ga]† = g, we h i have the following relation between the Bloch matrices for j <0. We now artificially extend the scattering regionby in- B†F(±)=Fa†(±)B. (42) cluding slices from the left and right leads and let it run fromj <0 to S+1+i with i>0. Eq.(30) gives for the transmission matrix element B. Equivalence of mode matching and Green function approaches τ′ =u† (+)G (E) G(0)(E) −1u (+) n,m R,n S+1+i,j j,j L,m h i (49) The retarded surface Green functions of the left and Using (47)eand (48) then gives right leads can be derived from Eqs. (40) and (41) τ′ = u† (+)Fi (+)G (E) g (E)=F−1(−) B† −1 ; g (E)=F (+)B−1. (43) n,m R,n R S+1,0 L L L R R R −1 G(0)(E) F−j(+)u (+) h i e 0,0 L L,m The retarded self-energies of Eq. (26) are then given by = λhi (+)λi−j (+)τ R,n L,m n,m Σ (E)=B F−1(−); Σ (E)=B†F (+). (44) = eiατ , (50) L L L R R R n,m Comparing Eqs. (13) and (22) then establishes withαreal. ThethirdlinefollowsfromapplyingEq.(7). The last line follows from the fact that we are only in- G′ (E)=G (E). (45) terestedinpropagatingstatesandforpropagatingstates i,j i,j 7 |λ| = 1. Using this result in (16) and (17) proves the FromEqs.(6)and(9)itistheneasytoshowthatthedual invariance of the total transmission probability with re- vectorsu (±)formthecolumnsofthematrix U(±)−1 † n spect to moving the boundaries between leads and scat- and that the F(±) matrices obey the equation tering region into the leads. (cid:2) (cid:3) e F(±)U(±)=U(±)Λ(±). (57) D. The Caroli expression In a similar way matrices Ua(±) and Λa(±) can be con- structed, see Eq. (36). ThetotaltransmissionprobabilityisgivenbyEq.(17), Usingthesedefinitions,onecangeneralizetheτ-matrix wherethesumhastobeoverpropagatingstatesonly. We of Eq. (20) to canextendthesummationtoincludeallN states(propa- gatingandevanescent)bydefininganN×N transmission τ =U−1(+)GrQ U (+). (58) R 0 L matrix Notethatτ-matrixelementsaredefinednotonlybetween 1 1 t=V2(+)τ V2(+), (51) propagating states, but also between evanescent states. R L However, as we remarked above already, only propagat- where τ is the matrix whoseeelements are given by ing states contribute to the physical transmission. Eq.(20). V (+)is defined asthe singular,diagonalma- R The second step is to express the velocity matrices in trix that has the velocities vR,n times the constant ~/aR termsofthe Γ-matrices. Todothis weuse anexpression on the diagonal for the right-propagating states and ze- for the velocity matrix, ros for evanescent states. We call it the velocity ma- trix. Likewise a pseudo-inverse velocity matrix VL(+) V(±)=i U†(±)B†U(±)Λ(±)−Λ†(±)U†(±)BU(±) , can be defined, which has 1/vL,n×aL/~ on the diagonal (59) for left-propagating states and all other matrix eleements which can(cid:2)be shown (see Appendix A for a proof) to b(cid:3)e are zero. These velocity matrices project onto the space equivalent to the definition introduced in the first para- of the propagating states so that the transmission ma- graph of this section. Using (57), this can be rewritten trixhasonlynon-zerovaluesbetweenpropagatingstates. for the right lead as (17) can then be expressed in the familiar form V (+) = iU†(+) B†F (+)−F†(+)B U (+) T = Tr t†t (52) R R R R R R R h i = Tr(cid:2)τ†V(cid:3)R(+)τVL(+) . = iU†R(+) ΣR−Σ†R UR(+) h i = U†(+)Γh U (+). i (60) Using the above definition of theevelocity matrix and R R R Ando’sexpressionsforthetransmissionmatrixelements, The second line follows from (44). A similar relation we will show how (52) can be rewritten as between the Γ-matrix and the velocity matrix for the left lead can be shown to exist, T =Tr[Γ GrΓ Ga], (53) R L where Gr,Ga are short-hand notations for GS+1,0(E) VL(+)=UaL†(−)ΓLUaL(−), (61) and Ga (E), respectively. The matrices Γ are de- 0,S+1 L/R fined as by using (42) and an equivalent expression for the ve- locity matrix, Eq. (A7). Eqs. (60) and (61) imply that Γ =i Σ −Σ† . (54) the Γ-matrices project onto the space spanned by the L/R L/R L/R propagating states. h i The third step is to introduce a matrix P that explic- Eq. (53) is known as the Caroli expression,15 and it is itly projects onto the propagating states of the left lead often used to calculate transmission probabilities.16,22 It is equivalent to the Kubo-Greenwood expression for the P = U (+)I [Ua(−)]−1 linear response regime.14,40 L p L The first step is to construct N ×N matrices U(±), Np the columns of which are the eigenmodes un(±), and = uL,n(+)uaL†,n(−), (62) diagonal matrices Λ(±), the elements of which are the n=1 X eigenvalues λn(±) e where the I -matrix contains 1 on the N diagonal ele- p p U(±)= u (±) u (±) ... u (±) , (55) mentsthatcorrespondtopropagatingstates,and0atall 1 2 N other positions. Given this projector matrix, it is possi- (cid:0) (cid:1) ble to prove that Λ(±) =λ (±)δ . (56) Q P=iΓ . (63) m,n n m,n 0 L 8 The proof is given in Appendix B. Using this property one has Q U (+)V (+)=Q PUa(−)V (+) 0 L L 0 L L −1 =iΓ Ua(−)V (+)=i Ua†(−) V (+)V (+) L L e L L e L L −1 h i =i Ua†(−) e I . e (64) L p h i FIG.2: Theparametersofaone-bandnearestneighbortight- Substituting (58), (60) and (64) into (52) leads directly bindingmodelforasingleimpurityinaonedimensionalchain. to the Caroli expression, Eq. (53). The roots λ(±) can be given a more familiar form. For E. Transmission matrix: a compact expression |ω|≤1 we define a wave number k by cos(ka)=ω, (70) Using the results of the previous section it is possible toderiveamorecompactexpressionforthetransmission where a is the lattice parameter. From Eqs. (69), (70) matrix. Combining (51), (58) and (62) one has one then obtains 1 1 λ(±)=e±ika, (71) t=V2(+)U−1(+)GrQ PUa(−)V2(+). (65) R R 0 L L which describes propagating states. For |ω| > 1 one de- By following the same steps as in Eq.e(64) and using fines κ by 1 1 V V2 =V2 this can be simplified to L L L cosh(κa)=|ω|. (72) 1 1 et=iVR2(+)U−R1(+)Gr[UaL†(−)]−1VL2(+). (66) One obtains λ(±) = exp(∓κa) if ω > 1 and λ(±) = −exp(∓κa) if ω < −1; both cases describe evanescent Writingoutthetransmissionmatrixelementsgivesthe states. compact expression Since the scattering region consists of a single im- purity, S = 1, Eqs. (12)-(14) give three linear equa- v v t =i~ R,n L,m u† (+)G (E)ua (−). (67) tionswiththreeunknownsdescribingthescatteringprob- n,m r aRaL R,n S+1,0 L,m lem. In a one-channel model one has F(±) = λ(±) and λ(±)−1 = λ(∓). There is only one possible incoming This is the tight-bindieng equivalent of tehe Fisher-Lee wave, so c (+) = 1. The linear equations then become expression relating transmission and Green function 0 in matrix form matrices.37 c β[λ(−)−λ(+)] 0 A c = 0 , (73) 1     VI. EXAMPLES c 0 2     with A. A simple analytical model E−ǫ−βλ(+) −β 0 L Weconsiderasystemconsistingofasingleimpurityin A= −β E−ǫ −β . L I R   aone-dimensionalchainandtreatthiswithinaone-band 0 −β E−ǫ−βλ(+) R nearestneighbortight-bindingmodel. Theparametersof  (74) this model are givenin Fig. 2. This model can be solved SolvingthissetofequationsandusingEqs.(70)and(71) analytically,23 soitcanserveasasimpletesttoillustrate we obtain the compact expression the equivalence of the different approaches. −ifsin(ka) We will first solvethe problem using the mode match- c =e2ika , (75) 2 ingapproachofSec.III. Itisconvenienttodefineascaled d+(1−b)cos(ka)−ibsin(ka) energy by defining the dimensionless parameters ω ≡ E2−βǫ. (68) b= βL2 +βR2, d= ǫ−ǫI, f = βLβR. (76) 2β2 2β β2 Themodelinvolvesonlyonechannelandwithum(±)=1 Applying Eqs. (15)-(17) then yields for the total trans- Eq. (4) reduces to mission probability −β+(E−ǫ)λ(±)−βλ(±)2 =0, (69) T(E)=|c |2. (77) 2 9 UsingEqs.(70)and(75)itiseasytoshowthatthistrans- setting k = iκ and setting the denominator to zero in mission probability is identical to Eq. (15) of Ref. 23, Eq. (75). This leads to the equation which was obtained using a different technique. It is instructive to solve the same problem using the (ω+d) ω+sgn(ω) ω2−1 −b=0; |ω|>1, (84) Green function approach of Sec. IVA. First one has (cid:16) p (cid:17) to find the surface Green functions of the leads from therootsofwhichgivetheenergiesofthelocalizedstates. Eqs.(24) and(25), which for the currentmodel become Againtheseresultsareequivalenttotheresultsobtained using the approach of Ref. 23. E−ǫ−β2g(E) g(E)=1, (78) WithintheGreenfunctionapproachtheenergiesofthe (cid:2) (cid:3) localizedstates aregivenby the poles ofthe Greenfunc- whereE isarealenergy. Thisequationhasthesolutions tion matrix. Via Eq. (81) this again leads to Eq. (84). Alternatively, since the Green function matrix is the in- e±ika verse of the A-matrix, cf. Eq. (80), its poles are given g±(E)= , (79) by the roots of det(A) = 0. This equation is equiva- β lent to Eq. (84), as is easily shown by setting λ(+) = ±exp(−κa) in the A-matrix. for both leads. The Green function matrix in the scat- tering regioncanthen be found fromEqs.(21)and (22), which can be combined in the 3×3 matrix equation B. Fe|vacuum|Fe tunnel junction AG(E)=I, (80) As anexample of a more complex system, we consider where G (E), i,j =0,...,2 are the matrix elements of an Fe|vacuum|Fe tunnel junction where the electronic i,j G(E)andAisgivenby Eq.(74). InvertingAyieldsthe structure is treated using the local density approxima- matrix element tion of DFT.42 The calculations are based upon a tight- bindingmuffintinorbital(TB-MTO)atomicspheresap- f e2ika proximation(ASA)implementation8,19,20,21,25 ofthefor- G (E)= , (81) 2,0 2βd+(1−b)cos(ka)−ibsin(ka) malism described in Sect. III. The first step in the calculation is the self-consistent with the parameters b,d and f defined by Eq.(76). The determination of the electronic structure of the tun- relevantGreenfunctionmatrixelementfortheideallead nel junction using the layer Green function approach of is found from Eq. (34) Ref. 35. The Fe leads are oriented in the (001) direc- tion and the atoms at the Fe(001) surfaces are kept at i G(0)(E)= . (82) their unrelaxed bulk positions. For the bcc structure 0,0 2βsin(ka) and TB-MTOs,43 a principal layerin the (001) direction contains two monolayers of Fe with a thickness of 2.866 Using these results in Eq. (30) one observes that the ex- ˚A. The vacuum region is modeled by a number of such pression for the (one channel) transmission matrix ele- slices, of the same thickness, filled with “empty” spheres ment becomes identical to Eq. (75). of the same size and packing as the Fe atomic spheres. Finally one can calculate the transmission probability The atomic sphere potentials of the vacuum region and fromtheCaroliexpressiongiveninSec.VD,cf. Eq.(53). four monolayers (two principal layers) of Fe on either Using Eqs. (26), (54) and (79) one obtains side of the vacuum are calculated self-consistently while the potentials of more distant layers are kept at their ΓL =ΓR =−2βsin(ka), (83) bulk values. These potentials then form the input to a transmission calculation based on mode matching.8,20,21 for left and right leads. Using Eqs. (81), (83) and Further technical details can be found in Ref. 25. Ga0,2 = (G2,0)∗ in Eq. (53) then yields an expression for A useful quantity to extract from the self-consistent the transmissionprobabilitythatisidenticaltoEq.(77). layer calculation is the layer density of states (LDOS) It illustrates the equivalence of the different approaches ρ (E). ItisrelatedtotheretardedGreenfunctionmatrix i for calculating the transmission in this simple model. defined in Eq. (21) by In addition to providing a channel for propagating states, an impurity can also give rise to localized states, ρ (E+iη)=−π−1ImTr[G (E+iη)] (85) i i,i whose energyis outside the energyband ofthe chain, cf. Eq. (72). Such a state does not contribute to the phys- wherethetracereferstotheusuallmangularmomentum ical transmission, but the transmission amplitude has a indices characterizing MTOs. For reasons of numerical poleatanenergythatcorrespondstoalocalizedstate.41 stabilitytheretardedGreenfunctionmatrixiscalculated Within the mode matchingapproachthis correspondsto by adding a finite imaginary partto the energy; we have anenergyatwhichc becomesinfinite. Forthepresent used η = 0.0025 Ry. In the mode matching approach, S+1 modelthe energiesoflocalizedstates canbe obtainedby theLDOScanbe directlyexpressedintermsofthe wave 10 FIG.3: (Color online) Layerdensityofstates (LDOS)at the Fe(001) surface layer. Top panel: majority spin. Top curve (red, dashed): as calculated using the layer Green function method.35 Bottomcurve(blue,solid): ascalculatedusingthe mode matching approach. Both calculations use η = 0.0025 Ry. For clarity, the top curve has been displaced by 10 units along the y-axis. Bottom panel: minority spin, the LDOS is shown with a negative sign. functionsor,alternatively,intermsoftheGreenfunction matrix of Eq. (18) ρ (E)=−π−1ImTr[G′ (E)]. (86) i i,i This Green function matrix can be calculated for a real energy, but in order to make a comparisonto the results obtained with Eq. (85), we add an imaginary part, η = 0.0025 Ry. In Fig. 3 we compare the LDOS obtained using (85) and (86) for the topmost Fe monolayer of (001) Fe|vacuum|Fe where the vacuum layerwasso thick (four principal layers, corresponding to a thickness of 11.466 ˚A) that the LDOS corresponds closely to that of a free Fe(001) surface. The two curves, displaced vertically for clarity, are indistinguishable, illustrating the equiv- alence of the two Green functions defined in (18) and (21). Moreover, the LDOS are in essential agreement with results found previously for a Fe(001) surface.44 If one resolves the LDOS into contributions from dif- ferent parts of the surface Brillouin zone, then the con- tribution at Γ (k = 0) exhibits sharp peaks near the k Fermi level. These are associated with characteristic FIG. 4: (Color online) The top panels represent the band surface states found on (001) surfaces of bcc transition structures along Γ-H of bulk Fe for majority and minority metals.45 These surface states are mainly derived from spins. Panels (a-h) give the ∆1 contribution to the LDOS d3z2−r2 orbitals on the surface atoms projecting into the ρi(E) at Γ of the Fe|vacuum interface as a function of the vacuum, and they have ∆ symmetry.46,47 They become layerposition. (a)LDOSofanFelayerataninfinitedistance 1 from the interface; the double peak structure below -4 eV most clearly visible if one filters out the contribution to and single peak (of a double peak structure) above -1 eV the LDOS at Γ of the states with ∆ symmetry. 1 correspond to the ∆1 bulk bands. (b) LDOS of the 6th Fe This contribution, calculated using the mode match- layer (the surface layer is layer 1). (c) LDOS of the 3rd Fe ing approach, is shown in Fig. 4 as a function of the layer; the sharp peak in the bulk band gap corresponds to distance from the surface. In the majorityspin LDOS, a a surface state. (d) LDOS of the the surface Fe layer; the sharppeak isfoundatE −1.8eV,whichisverypromi- surface state peak has maximum amplitude. (e-h) Same for F nent in the surface Fe layer, cf. Fig. 4(d) and decays minority spin.

ebook img

Conductance calculations for quantum wires and interfaces: mode matching and Green functions PDF

0.38 MB·English
by  P. A. Khomyakov
#additional_collections #journals #arxiv
Save to my drive
Quick download
Download
Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview Conductance calculations for quantum wires and interfaces: mode matching and Green functions

Conductance calculations for quantum wires and interfaces: mode matching and Green functions P. A. Khomyakov, G. Brocks,∗ V. Karpan, M. Zwierzycki, and P. J. Kelly Computational Materials Science, Faculty of Science and Technology and MESA+ Research Institute, University of Twente, P. O. Box 217, 7500 AE Enschede, The Netherlands. 5 (Dated: February 2, 2008) 0 Landauer’s formularelates theconductanceofaquantumwire orinterfacetotransmission prob- 0 abilities. Totaltransmission probabilitiesarefrequentlycalculated usingGreenfunctiontechniques 2 and an expression first derived by Caroli. Alternatively, partial transmission probabilities can be n calculated from the scattering wave functions that are obtained by matching thewave functions in a the scattering region to the Bloch modes of ideal bulk leads. An elegant technique for doing this, J formulated originally by Ando, is here generalized to any Hamiltonian that can be represented in 5 tight-bindingform. Amorecompactexpressionforthetransmissionmatrixelementsisderivedand 2 it is shown how all the Green function results can be derived from the mode matching technique. We illustrate this for a simple model which can be studied analytically, and for an Fe|vacuum|Fe ] tunneljunction which we study using first-principlescalculations. i c s PACSnumbers: 73.63.-b,73.40.-c,73.20.-r,85.35.-p - l r t m I. INTRODUCTION matching17 and Green function15,16,22 approaches is not immediately obvious. Indeed, it wasrecently statedthat . at Since the discovery of giant magnetoresistance in Ando’s approach is incomplete and does not yield the m metallic multilayers there has been considerable interest correct expression for the conductance.22 - in studying electronic transport in layered materials.1,2 d At the same time, experimental control of the lateral n scalehasenabledstudies ofelectronictransportinquan- Inthis paper we demonstrate that the twoapproaches o tum wires of atomic dimensions.3 Because of the small are completely equivalent. In the Green function ap- c [ dimensions involved, the transport properties of both of proach,asmallimaginarypartmustbe addedtoorsub- these systems should be understoodon the basis oftheir tracted from the energy in order to distinguish between 1 atomic structure. This perception has generated a large the retarded and advanced forms.5,6,7,10,11,12,13,14,22 In v effortinrecentyearstocalculatetheconductanceofmul- mode matching,scatteringwavefunctions arecalculated 9 0 tilayers and quantum wires from first principles. Several thatincorporatetheretardedoradvancedboundarycon- 6 different approaches have been formulated which have a ditions directly. This makes it possible to solve the 1 common basis in the Landauer-Bu¨ttiker approachor are scattering problem also at real, instead of complex en- 0 equivalent to it. In the linear response regime, the con- ergies. In addition to yielding the total conductance, by 5 ductance G is expressed as a quantum mechanical scat- focussingonwavefunctionsthecontributionofeachindi- 0 tering problem4 and can be simply related to the total vidualscatteringchannelcanbeidentified. Inparticular, / t transmissionprobabilityatthe Fermienergy,T(E ), as we derive a simple, compactexpressionfor the transmis- a F m sion matrix elements, see Eq. (67). e2 G = T(E ). (1) - F h d n The multilayer or quantum wire is generally considered The paper is organized as follows. In the next section o as a scattering region of finite size, sandwiched between the Hamiltonian we will use is introduced. This model c two semi-infinite ballistic wires. Aiming at a materials allows us to study both quantum wires that are finite : v specific description, most current approaches treat the in the directions perpendicular to the wire, and systems i X electronicstructurewithintheframeworkofdensityfunc- that are periodic in these directions such as single inter- tional theory (DFT).5,6,7,8,9,10,11,12,13,14 faces,sandwichesandmultilayers. We willuse the single r a FrequentlytheconductanceiscalculatedusingaGreen term “quantum wire” to describe both systems. In Secs. function expression first derived by Caroli et al.15,16 III and IV the mode matching and Green function tech- An alternative technique, suitable for Hamiltonians that niques are summarized. The equivalence of the trans- can be represented in tight-binding form, has been for- mission matrices obtained using these two techniques is mulated by Ando.17 It is based upon directly match- demonstratedinSec. VandtheCaroliexpressionforthe ing the wave function in the scattering region to the conductance is derived from the mode matching expres- Bloch modes of the leads. The latter technique has sions. InSec. VIthetwotechniquesareappliedfirsttoa been applied to conductance calculations at the empiri- simple analytical model,23 and then to an Fe|vacuum|Fe caltight-binding level,18 as wellas onthe first-principles tunnel junction using numerical first-principles calcula- DFT level.8,19,20,21 The relationship between the mode tions. ThemainconclusionsaresummarizedinSec. VII. 2 II. HAMILTONIAN Wesetupatight-bindingrepresentationofthe Hamil- tonian. This is not a severe restriction since a first- principles DFT implementation that uses a representa- tion on an atomic orbital basis set has the same math- ematical structure as a tight-binding model.24,25 Alter- natively, an implementation that uses a representation of the Hamiltonian on a grid in real space, can also be mappedontoatight-bindingmodel.26Webeginbydivid- FIG. 1: (Color online) Hamilton matrix of a quantum wire divided into slices. The left (L) and right (R) leads are ideal ing the system into slices (“principal layers”)perpendic- periodic wires that span the cells i = −∞,...,0 and i = ular to the wire direction.27 The thickness of these slices S+1,...,∞, respectively. The scattering region spans cells is chosen such that there is only an interaction between i=1,...,S. neighboring slices. Labelling each slice with an index i, the Schr¨odinger equation of the quantum wire becomes A. Bloch matrices −H c +(EI−H )c −H c =0, (2) i,i−1 i−1 i,i i i,i+1 i+1 Thefirststepconsistsoffindingsolutionsfortheleads, fori=−∞,...,∞. AssumingthateachslicecontainsN for which Eq. (2) can be simplified to sites and/or orbitals, c is a vector of dimension N con- i taining the wave function coefficients on all sites and/or −Bc +(EI−H)c −B†c =0. (3) i−1 i i+1 orbitals of slice i. The N ×N matrices H and H i,i i,i±1 consist of on-slice and hopping matrix elements of the These equations hold for i = −∞,...,−1 and i = Hamiltonian, respectively. I is the N ×N identity ma- S +2,...,∞, i.e. the left and right leads. To keep the trix. A schematic representation of the structure of the notation as simple as possible, we have omitted the sub- Hamiltonian is given in Fig. 1. scripts L and R, see Fig. 1. Since the leads are peri- Eq. (2) is valid both for quantum wires that are finite odic wires, one can make the ansatz that the solutions in the directions parallel to the slices, and for layered have Bloch symmetry, i.e. ci = λci−1, ci+1 = λ2ci−1, systems that are periodic in these directions. In the lat- whereλistheBlochfactor. SubstitutingthisintoEq.(3) ter case, translations in the transverse direction can be results in a quadratic eigenvalue equation of dimension described in terms of a Bloch wave vector in the two- N. The latter canbe solvedmost easilyby transforming dimensionalBrillouinzone,k ,whichisagoodquantum it to an equivalent linear (generalized) eigenvalue prob- k number and the system becomes effectively one dimen- lem of dimension 2N and solving this by a standard sional. Explicit expressions for the Hamilton matrix el- algorithm.29,30 ements depend upon the particular localized orbital ba- It can be shown that the equation generally has 2N sis or realspace grid representationused.28 Since details solutions,whichcanbedividedintoN right-goingmodes of the tight-binding muffin tin orbital scheme used in and N left-going modes,31 labelled as (+) and (−) in Refs. 8,19,20,21 are given in Ref. 25 and of the real- the following. Right-going modes are either evanescent space high-order finite difference method can be found wavesthataredecayingtotheright,orwavesofconstant in Ref. 26, they will not be discussed further here. amplitudethatarepropagatingtotheright,whereasleft- The system is divided into three parts, with i = goingmodes aredecaying or propagatingto the left. We −∞,...,0correspondingtotheleftlead(L),i=1,...,S denote the eigenvalues by λn(±) where n = 1,...,N, to the scattering region (S) and i = S + 1,...,∞ to the corresponding eigenvectors by un(±) and write the the right lead (R). The leads are assumed to be ideal eigenvalue equation as wires characterized by a periodic potential. It is then −Bu (±)+(EI−H)λ (±)u (±)−B†λ (±)2u (±)=0. n n n n n sufficient to identify a slice with a translational period (4) along the wire. By construction, the Hamilton matrix is In the following we assume that the vectors u (±) are the same for each slice of the leads, i.e. H ≡ H , n i,i L/R normalized;notethatingeneraltheyarenot orthogonal. Hi,i−1 ≡ BL/R and Hi,i+1 ≡ B†L/R for the left/right One can distinguish right- from left-going evanescent leads. Fig. 1 summarizes our model of a quantum wire. modes on the basis of their eigenvalues; right going evanescentmodes have |λ(+)|<1 and left going evanes- cent modes have |λ(−)| > 1. Propagating modes have III. MODE MATCHING APPROACH |λ(±)| = 1, so here one has to determine the Bloch ve- locityanduseitssigntodistinguishrightfromleftprop- Eq. (2) can be solvedin a particularly convenient way agation. TheBlochvelocitiesaregivenbytheexpression by a technique which we call mode matching. In this, we follow Ando and generalize his approach to a slice 2a geometry.17 vn(±)=− ~ Im λn(±)un(±)†B†un(±) , (5) (cid:2) (cid:3) 3 where a is the slice thickness. A derivation of this equa- with a modified Hamilton matrix defined so H′ =H i,j i,j tion is given in Appendix A. for the elements {i,j =0,1}, {i=1,...,S; j =i,i±1} Since the eigenvectors are non-orthogonal,it is conve- and {i,j =S+1,S}, but with nient to define dual vectors u (±) n H′ = H +B F−1(−), u†(±)u (±)=δ ; u†(±)u (±)=δ . (6) 0,0 L L L n m n,m e n m n,m H′ = H +B†F (+). (13) S+1,S+1 R R R Any wave function in the leads can be expressed as a e e linear combination of the lead modes. This can be done H′ =0 andH′ =0,so the modifiedscattering 0,−1 S+1,S+2 in a compact form using the two N ×N Bloch matrices regionisdecoupledfromtheleads. Ontheleft-handside for right- and left-going solutions of Eq. (12), we have a source term with N Q =B F−1(+)−F−1(−) , (14) F(±)= λ (±)u (±)u†(±). (7) 0 L L L n n n n=1 (cid:2) (cid:3) X and Q = 0 for i = 1,...,S +1. Eq. (12) defines a set i For any integer i, Fi is givenby a simeilar expressionbut of (S +2)×N linear equations. Because the Hamilton with λ replaced by λi. From the foregoing it is easy to matrixisblocktridiagonal,eachblockbeingofdimension n n show that the Bloch matrices obey the equation N, this set of equations can be solved efficiently using a block Gaussian elimination scheme.30 The total wave −BF−1(±)+(EI−H)−B†F(±)=0. (8) function c can then be obtained by back substitution. i The transmission is obtained from the wave function A general solution in the leads can be expressed in in the right lead c (+). In particular, choosing the terms of a recursion relation S+1 incoming wave as one of the modes of the left lead, i.e. ci = ci(+)+ ci(−) c0(+) = uL,m(+), generalized transmission matrix ele- = Fi−j(+)c (+)+Fi−j(−)c (−). (9) ments τn,m are defined by expanding cS+1(+) in modes j j of the right lead Fixingthecoefficientsinoneslicethensetstheboundary conditionsanddeterminesthesolutionforthewholelead. N c (+)= u (+)τ . (15) S+1 R,n n,m n=1 X B. Transmission matrix Byletting c (+) runoverallpossible incomingmodes of 0 the left lead u (+); m = 1,...,N, a full transmission The scattering region is defined by i = 1,...,S, see L,m matrix is obtained. Fig. 1. Right and left of the scattering regions one has Matrix elements can be defined for all modes, prop- the recursion relations for the states in the leads from agating and evanescent, but of course only matrix ele- Eq. (9) ments where n,m denote propagating modes contribute c = F−1(+)c (+)+F−1(−)c (−) to the real physical transmission. These modes can be −1 L 0 L 0 = F−1(+)−F−1(−) c (+)+F−1(−)c ,(10) selected by making use of their eigenvalues; see the dis- L L 0 L 0 cussionfollowingEq.(4). Thephysicaltransmissionma- with c = c(cid:2)(+)+ c (−) and (cid:3) trixelementsarethenfoundbynormalizingwithrespect 0 0 0 to the current16 c =F (+)c (+)+F (−)c (−), (11) S+2 R S+1 R S+1 v (+)a where the subscripts L and R distinguish between the t = R,n L τ , (16) n,m n,m Bloch matrices of the left and right leads. svL,m(+)aR The boundary conditions for the scattering problem are set up in the usual way. The vector c0(+) is treated where vL,m(+) and vR,n(+) are the Bloch velocities in as the source, i.e. as a fixed incoming wave from the left the directionofthe wireforthe right-propagatingmodes lead. There is no incoming wave from the right lead, so m and n in the left and right leads, respectively, see we set cS+1(−)=0. Eq. (5); aL and aR are the slice thicknesses of left and Havingsettheboundaryconditions,Eqs.(10)and(11) right leads.32 The total transmissionprobability is given can be used to rewrite Eq. (2) in the region not covered by by Eq. (3), i.e. for i = 0,...,S +1. This describes the wave function in the scattering region and the matching (+) to the solutions in the leads. Eq. (2) in this region is T(E)= |tn,m|2, (17) rewritten as n,m X −H′ c + EI−H′ c −H′ c =Q c (+), and the conductance is given by Eq. (1) evaluated at i,i−1 i−1 i,i i i,i+1 i+1 i 0 (12) E =E . F (cid:0) (cid:1) 4 C. Green function matrix Here g (z) and g (z) are the surface Green functions of L R the semi-infinite left and right leads, respectively, which can be calculated using an iterative technique: denoting Solving the set of linear equations Eq. (12) directly leads to the conductance. However, to facilitate a con- G[n](z)asthesolutionofanequationsimilartoEq.(21), i,j nection to the Green function approach discussed in but with H = 0 for {i > n∨j > n}, one can easily i,j the next section, we can formulate the solution in a derive the right-going recursion relation slightly different way. A finite Green function matrix G′ (z), i,j =0,...,S+1 can be defined by zI−Hn+1,n+1−Hn+1,nG[nn,]n(z)Hn,n+1 i,j h G[n+1] (z)=Ii. (23) −H′ G′ + zI−H′ G′ −H′ G′ =Iδ , n+1,n+1 i,i−1 i−1,j i,i i,j i,i+1 i+1,j i,j (18) Foranidealwirewithi,j =−∞,...,n,G[n] (z)=g (z) with z complex. (cid:0)Note, how(cid:1)ever, that the matrices H′0,0 should be independent of n resulting inn,tnhe followLing and H′S+1,S+1 are non-Hermitian and G′i,j(E) is also equation for the surface Green function, uniquely defined for real energies. The Green function matrix allows the solution of Eq. (12) to be written as zI−HL−BLgL(z)B†L gL(z)=I. (24) h i c =G′ (E)Q c (+). (19) Several iterative algorithms have been formulated for i i,0 0 0 solving this non-linear matrix equation.27,35,36 A simi- Asbefore,thetransmissioncanbeextractedati=S+1 lar reasoning based upon a left-going recursion for the and comparison with Eq. (15) gives right lead results in an equation for the surface Green function g (z) of the right lead R τn,m =u†R,n(+)G′S+1,0(E)Q0uL,m(+). (20) zI−H −B†g (z)B g (z)=I. (25) R R R R R which can be useed in Eq. (16). The Green function ma- Again, sethting z = E +iη in Eqs.i(24) and (25) defines trix block G′ (E) can be found by solving Eq. (18) the usualretardedsurface Greenfunctions g (E). Al- S+1,0 L/R using a recursive algorithm that resembles a Gaussian though we are mainly interested in the physical limit elimination scheme.25,33 lim , in practice a finite value of η is often used in η→0 order to make the iterative algorithms stable. The quantities IV. GREEN FUNCTION APPROACH Σ (E)=B g (E)B† ; Σ (E)=B†g (E)B , (26) L L L L R R R R An apparently quite different route to the transmis- which appear in Eqs. (22)-(25), are called the self- sion matrix starts by defining an infinite Green function energies of the left and right leads, respectively.16 Once matrix G (z) for i,j =−∞,...,∞ with respect to the theseareobtained,thefiniteHamiltonmatrixofEq.(22) i,j original Hamiltonian of Eq. (2). is constructed and the retarded Green function matrix G (E) can be found using a recursive algorithm.33 Us- i,j −H G +(zI−H )G −H G =Iδ . ingtheleadmodesthe transmissionmatrixelementscan i,i−1 i−1,j i,i i,j i,i+1 i+1,j i,j (21) then be calculated, as will be shown in the next section, Choosing z = lim (E ±iη) defines as usual the re- Sec.IVB. Alternatively,thetotaltransmissionprobabil- η→0 tarded/advanced Green function matrix. We shall use ity can be expressed in a form that does not require the G (E) to denote the retarded Green function matrix lead modes explicitly, which is discussed in Sec. V. i,j and Ga (E) to denote the advanced Green function ma- i,j trix. B. Transmission matrix A. Partitioning The transmission matrix can be obtained from the Green function matrix of Eq. (21). To do this, we adaptaFisher-Leetypeofapproachtoourtight-binding Eq.(21)ismostconvenientlysolvedbyapplyingapar- formulation.37 Assuming that the unperturbed reference titioningtechnique.16,34Itisstraightforwardtoshowthat wave function is the Bloch mode u (+) that comes in thefinitepartG (z)definedfori,j =0,...,S+1canbe L,m i,j from the left lead, the Lippmann-Schwinger equation38 derivedfromaclosedsetofequations,similartoEq.(21), in tight-binding form is butwithH replacedbyH′′ whereH′′ =H forthe i,j i,j i,j i,j elements {i,j = 0,1}, {i = 1,...,S; j = i,i±1} and c = u (+)+ G V u (+) i L,m,i i,j j,k L,m,k {i,j =S+1,S}, but with j,k X H′′ (z) = H +B g (z)B†, 0,0 L L L L = Fi(+)+ G V Fk(+) u (+).(27) H′S′+1,S+1(z) = HR+B†RgR(z)BR. (22)  L Xj,k i,j j,k L  L,m   5 Hereu (+)isthereferencewavefunctioninslicei. It terms of its eigenmodes. The columns of such a Green L,m,i obeys Bloch symmetry and u (+) ≡ u (+) is the function obey the equation L,m L,m,0 Blochmodeattheorigin,seeSec.IIIA. ThematrixV j,k represents the perturbationwith respect to the ideal left −BG(0) + E+I−H G(0)−B†G(0) =Iδ , (31) i−1,j i,j i+1,j i,j lead. Eq. (27) can be simplified using the Dyson equation, where E+ = E(cid:0) + iη. F(cid:1)or i 6= j the solution is simi- which in tight-binding form reads lar to that of the wave functions, see Eq. (3). In ad- dition, the retarded Green function should consist only G = G(0)+ G V G(0) i,0 i,0 i,j j,k k,0 of propagating waves that move outwards from the δ- Xj,k source and/or evanescent states that decay away from the source.38 From Eq. (9), we have the column recur- = Fi(+)+ G V Fk(+) G(0), (28) sion relations  L i,j j,k L  0,0 j,k X G(0)(E) = Fi−j(−)G(0)(E), i<j,   i,j j,j using Eq. (32). Comparing Eqs. (27) and (28) one finds the simple expression G(i,0j)(E) = Fi−j(+)G(j0,j)(E), i>j. (32) −1 ci =Gi,0(E) G(00,0)(E) uL,m(+). (29) ThediagonalblockG(j0,j)(E)cannowbeobtainedbycom- h i bining Eqs. (31) and (32), which gives for i=j Fromthedefinitionofthegeneralizedtransmissionma- trix elements, cf. Eq. (15), one then obtains the expres- −BF−1(−)+E+I−H−B†F(+) G(0) =I. (33) sion j,j τ =u† (+)G (E) G(0)(E) −1u (+). (30) Elim(cid:2)inating E+I−H using Eq. (8) th(cid:3)en yields n,m R,n S+1,0 0,0 L,m h i −1 To find τne,m one needs to calculate only the Green func- G(j0,j)(E) =B F−1(+)−F−1(−) , (34) tion matrix blocks G (E) of the full system and S+1,0 h i h i G(0)(E) of the idealleft lead. The physicaltransmission or the equivalent 0,0 matrix elements and the total transmission probability can then be obtained from Eqs. (16) and (17). G(0)(E) −1 =B†[F(−)−F(+)], (35) j,j h i Eqs. (32) and (34) represent the full expression for the V. MODE MATCHING VERSUS GREEN FUNCTIONS Green function G(0) of an infinite ideal wire in terms i,j of the Bloch matrices F(±) and thus in terms of the The two seemingly different formalisms introduced in eigenmodes. Note that we can set E+ = E since, in Secs. IIIandIVareinfactcloselyrelated. Inthissection terms of the modes, the retarded Green function matrix we will show how all Green function results can be ob- is uniquely defined for real energies. tained from the mode matching approach. We begin by The advanced Green function matrix G(0)a(E) can be i,0 expressing the Green function matrices of ideal wires in found from a similar procedure. It should consist of termsoftheBlochmatrices,F(±). Theseexpressionsare propagatingwaves that move towards the source and/or thenusedtoprovethatthetransmissionmatrixelements evanescentstatesthatgrowtowardsthesource. Onecan obtainedfromthemodematchingandGreenfunctionap- constructtwonewBlochmatricesFa(±),whicharesim- proaches,cf. Eqs.(20)and(30),areinfactidentical. Af- ilar to those defined in Eq. (7). In Fa(+) one collects ter that, we show that the transmission matrix elements the modes thatare decayingto the right(growingto the are independent of the exact positions within the leads left) and modes that are propagating to the left. Fa(−) thatareusedtomatchtheleadstothescatteringregion, thencontainsmodesthatgroworpropagatetotheright. apart from a trivial phase factor. Then we derive from the mode matching expressionfor the total transmission N Fa(±) = λa(±)ua(±)ua†(±), with (36) probability the Green function expression known as the n n n Caroli expression.15 Finally, a more compact expression nX=1 for the transmission matrix elements is derived. λan(±) = λn(∓), uan(±)=eun(∓)propagating λa(±) = λ (±), ua(±)=u (±)evanescent. n n n n A. Green functions of ideal wires in terms of Bloch Using these definitions, expressions for the advanced matrices Green function matrix are obtained from Eqs. (32) and (34) by replacing F(±) with Fa(±). We begin by deriving an expression for the retarded From the general relation between retarded and ad- Green function matrix G(0) of an ideal infinite wire in vanced Green functions, G = Ga †, the following i,j i,j j,i (cid:0) (cid:1) 6 row recursion relations can be deduced for the retarded Inotherwords,thetwoGreenfunctionsdiscussedinSecs. Green function III and IV are identical. By comparing Eqs. (14) and (34) one has G(0)(E) = G(0)(E) Fa†(+) j−i, i<j, i,j i,i −1 G(0)(E) = G(0)(E)(cid:2)Fa†(−)(cid:3)j−i, i>j. (37) Q0 = G(00,0)(E) . (46) i,j i,i h i TheretardedGreenfunctio(cid:2)nG(s)(E(cid:3))ofasemi-infinite In conclusion, the two expressions for the generalized i,0 transmission matrix elements, Eqs. (20) and (30), are wire extendingfromi=−∞,...,0canbeobtainedusing identical. a similar technique. Instead of Eq. (32), we get G(s)(E)=Fi(−)g(E), i<0. (38) i,0 C. Invariance of transmission probability where g(E) = G(s)(E) is the surface Green function. 0,0 Apartfromtrivialphase factors,the transmissionma- Using this in Eq. (31) gives for i = 0 and j = 0 and for trix elements should not depend on where in the ideal i=−1 and j =0, respectively, lead the wave function matching is carried out. In a re- centpaperit wasstatedthat Ando’sexpressionfort , −BF−1(−)+E+I−H g=I, n,m Eqs. (20) and (16), lacks this invariance property and is (cid:2)−BF−1(−)+E+I−H(cid:3)F−1(−)g =B†g. (39) therefore incomplete.22 One can however prove directly from Eq. (20) that the transmission matrix elements do Note tha(cid:2)t the B† term is absen(cid:3)t in the first equation have the required invariance property.39 The proof be- since we aredealing with a semi-infinite wire. These two comes easier if the equivalence of Eqs. (20) and (30), equations can be easily solved to find an expression for established above, is used. the surface Green function The scattering region runs from slices 0 to S + 1 if we include the boundaries with the left and right leads. g(E)=F−1(−) B† −1. (40) This means that the Green function matrix G with i,j indicesi,j outsidethisregionobeystheequationsforthe Eqs.(38)and(40)representtheG(cid:0)ree(cid:1)nfunctionofasemi- idealleads. Fromthecolumnandrowrecursionrelations, infinite ideal wire extending from i = −∞,...,0. In Eqs. (32) and (37), one derives a similar fashion, one gets for the Green function of a semi-infinite ideal wire extending from i=0,...,∞ G (E)=Fi (+)G (E) Fa†(−) j, (47) S+1+i,j R S+1,0 L G(s)(E)=Fi+1(+)B−1, i≥0. (41) h i i,0 for j < 0, i > 0. In a similar way, one derives for the Green function matrix of the left lead Analogously to Eqs. (38)-(41), one can define the ad- vanced Green function matrix Gi(,s0)a(E) in terms of the G(0)(E) = Fj(+)G(0)(E) Fa†(−) j, (48) j,j L 0,0 L Bloch matrices Fa(±). Moreover, since [ga]† = g, we h i have the following relation between the Bloch matrices for j <0. We now artificially extend the scattering regionby in- B†F(±)=Fa†(±)B. (42) cluding slices from the left and right leads and let it run fromj <0 to S+1+i with i>0. Eq.(30) gives for the transmission matrix element B. Equivalence of mode matching and Green function approaches τ′ =u† (+)G (E) G(0)(E) −1u (+) n,m R,n S+1+i,j j,j L,m h i (49) The retarded surface Green functions of the left and Using (47)eand (48) then gives right leads can be derived from Eqs. (40) and (41) τ′ = u† (+)Fi (+)G (E) g (E)=F−1(−) B† −1 ; g (E)=F (+)B−1. (43) n,m R,n R S+1,0 L L L R R R −1 G(0)(E) F−j(+)u (+) h i e 0,0 L L,m The retarded self-energies of Eq. (26) are then given by = λhi (+)λi−j (+)τ R,n L,m n,m Σ (E)=B F−1(−); Σ (E)=B†F (+). (44) = eiατ , (50) L L L R R R n,m Comparing Eqs. (13) and (22) then establishes withαreal. ThethirdlinefollowsfromapplyingEq.(7). The last line follows from the fact that we are only in- G′ (E)=G (E). (45) terestedinpropagatingstatesandforpropagatingstates i,j i,j 7 |λ| = 1. Using this result in (16) and (17) proves the FromEqs.(6)and(9)itistheneasytoshowthatthedual invariance of the total transmission probability with re- vectorsu (±)formthecolumnsofthematrix U(±)−1 † n spect to moving the boundaries between leads and scat- and that the F(±) matrices obey the equation tering region into the leads. (cid:2) (cid:3) e F(±)U(±)=U(±)Λ(±). (57) D. The Caroli expression In a similar way matrices Ua(±) and Λa(±) can be con- structed, see Eq. (36). ThetotaltransmissionprobabilityisgivenbyEq.(17), Usingthesedefinitions,onecangeneralizetheτ-matrix wherethesumhastobeoverpropagatingstatesonly. We of Eq. (20) to canextendthesummationtoincludeallN states(propa- gatingandevanescent)bydefininganN×N transmission τ =U−1(+)GrQ U (+). (58) R 0 L matrix Notethatτ-matrixelementsaredefinednotonlybetween 1 1 t=V2(+)τ V2(+), (51) propagating states, but also between evanescent states. R L However, as we remarked above already, only propagat- where τ is the matrix whoseeelements are given by ing states contribute to the physical transmission. Eq.(20). V (+)is defined asthe singular,diagonalma- R The second step is to express the velocity matrices in trix that has the velocities vR,n times the constant ~/aR termsofthe Γ-matrices. Todothis weuse anexpression on the diagonal for the right-propagating states and ze- for the velocity matrix, ros for evanescent states. We call it the velocity ma- trix. Likewise a pseudo-inverse velocity matrix VL(+) V(±)=i U†(±)B†U(±)Λ(±)−Λ†(±)U†(±)BU(±) , can be defined, which has 1/vL,n×aL/~ on the diagonal (59) for left-propagating states and all other matrix eleements which can(cid:2)be shown (see Appendix A for a proof) to b(cid:3)e are zero. These velocity matrices project onto the space equivalent to the definition introduced in the first para- of the propagating states so that the transmission ma- graph of this section. Using (57), this can be rewritten trixhasonlynon-zerovaluesbetweenpropagatingstates. for the right lead as (17) can then be expressed in the familiar form V (+) = iU†(+) B†F (+)−F†(+)B U (+) T = Tr t†t (52) R R R R R R R h i = Tr(cid:2)τ†V(cid:3)R(+)τVL(+) . = iU†R(+) ΣR−Σ†R UR(+) h i = U†(+)Γh U (+). i (60) Using the above definition of theevelocity matrix and R R R Ando’sexpressionsforthetransmissionmatrixelements, The second line follows from (44). A similar relation we will show how (52) can be rewritten as between the Γ-matrix and the velocity matrix for the left lead can be shown to exist, T =Tr[Γ GrΓ Ga], (53) R L where Gr,Ga are short-hand notations for GS+1,0(E) VL(+)=UaL†(−)ΓLUaL(−), (61) and Ga (E), respectively. The matrices Γ are de- 0,S+1 L/R fined as by using (42) and an equivalent expression for the ve- locity matrix, Eq. (A7). Eqs. (60) and (61) imply that Γ =i Σ −Σ† . (54) the Γ-matrices project onto the space spanned by the L/R L/R L/R propagating states. h i The third step is to introduce a matrix P that explic- Eq. (53) is known as the Caroli expression,15 and it is itly projects onto the propagating states of the left lead often used to calculate transmission probabilities.16,22 It is equivalent to the Kubo-Greenwood expression for the P = U (+)I [Ua(−)]−1 linear response regime.14,40 L p L The first step is to construct N ×N matrices U(±), Np the columns of which are the eigenmodes un(±), and = uL,n(+)uaL†,n(−), (62) diagonal matrices Λ(±), the elements of which are the n=1 X eigenvalues λn(±) e where the I -matrix contains 1 on the N diagonal ele- p p U(±)= u (±) u (±) ... u (±) , (55) mentsthatcorrespondtopropagatingstates,and0atall 1 2 N other positions. Given this projector matrix, it is possi- (cid:0) (cid:1) ble to prove that Λ(±) =λ (±)δ . (56) Q P=iΓ . (63) m,n n m,n 0 L 8 The proof is given in Appendix B. Using this property one has Q U (+)V (+)=Q PUa(−)V (+) 0 L L 0 L L −1 =iΓ Ua(−)V (+)=i Ua†(−) V (+)V (+) L L e L L e L L −1 h i =i Ua†(−) e I . e (64) L p h i FIG.2: Theparametersofaone-bandnearestneighbortight- Substituting (58), (60) and (64) into (52) leads directly bindingmodelforasingleimpurityinaonedimensionalchain. to the Caroli expression, Eq. (53). The roots λ(±) can be given a more familiar form. For E. Transmission matrix: a compact expression |ω|≤1 we define a wave number k by cos(ka)=ω, (70) Using the results of the previous section it is possible toderiveamorecompactexpressionforthetransmission where a is the lattice parameter. From Eqs. (69), (70) matrix. Combining (51), (58) and (62) one has one then obtains 1 1 λ(±)=e±ika, (71) t=V2(+)U−1(+)GrQ PUa(−)V2(+). (65) R R 0 L L which describes propagating states. For |ω| > 1 one de- By following the same steps as in Eq.e(64) and using fines κ by 1 1 V V2 =V2 this can be simplified to L L L cosh(κa)=|ω|. (72) 1 1 et=iVR2(+)U−R1(+)Gr[UaL†(−)]−1VL2(+). (66) One obtains λ(±) = exp(∓κa) if ω > 1 and λ(±) = −exp(∓κa) if ω < −1; both cases describe evanescent Writingoutthetransmissionmatrixelementsgivesthe states. compact expression Since the scattering region consists of a single im- purity, S = 1, Eqs. (12)-(14) give three linear equa- v v t =i~ R,n L,m u† (+)G (E)ua (−). (67) tionswiththreeunknownsdescribingthescatteringprob- n,m r aRaL R,n S+1,0 L,m lem. In a one-channel model one has F(±) = λ(±) and λ(±)−1 = λ(∓). There is only one possible incoming This is the tight-bindieng equivalent of tehe Fisher-Lee wave, so c (+) = 1. The linear equations then become expression relating transmission and Green function 0 in matrix form matrices.37 c β[λ(−)−λ(+)] 0 A c = 0 , (73) 1     VI. EXAMPLES c 0 2     with A. A simple analytical model E−ǫ−βλ(+) −β 0 L Weconsiderasystemconsistingofasingleimpurityin A= −β E−ǫ −β . L I R   aone-dimensionalchainandtreatthiswithinaone-band 0 −β E−ǫ−βλ(+) R nearestneighbortight-bindingmodel. Theparametersof  (74) this model are givenin Fig. 2. This model can be solved SolvingthissetofequationsandusingEqs.(70)and(71) analytically,23 soitcanserveasasimpletesttoillustrate we obtain the compact expression the equivalence of the different approaches. −ifsin(ka) We will first solvethe problem using the mode match- c =e2ika , (75) 2 ingapproachofSec.III. Itisconvenienttodefineascaled d+(1−b)cos(ka)−ibsin(ka) energy by defining the dimensionless parameters ω ≡ E2−βǫ. (68) b= βL2 +βR2, d= ǫ−ǫI, f = βLβR. (76) 2β2 2β β2 Themodelinvolvesonlyonechannelandwithum(±)=1 Applying Eqs. (15)-(17) then yields for the total trans- Eq. (4) reduces to mission probability −β+(E−ǫ)λ(±)−βλ(±)2 =0, (69) T(E)=|c |2. (77) 2 9 UsingEqs.(70)and(75)itiseasytoshowthatthistrans- setting k = iκ and setting the denominator to zero in mission probability is identical to Eq. (15) of Ref. 23, Eq. (75). This leads to the equation which was obtained using a different technique. It is instructive to solve the same problem using the (ω+d) ω+sgn(ω) ω2−1 −b=0; |ω|>1, (84) Green function approach of Sec. IVA. First one has (cid:16) p (cid:17) to find the surface Green functions of the leads from therootsofwhichgivetheenergiesofthelocalizedstates. Eqs.(24) and(25), which for the currentmodel become Againtheseresultsareequivalenttotheresultsobtained using the approach of Ref. 23. E−ǫ−β2g(E) g(E)=1, (78) WithintheGreenfunctionapproachtheenergiesofthe (cid:2) (cid:3) localizedstates aregivenby the poles ofthe Greenfunc- whereE isarealenergy. Thisequationhasthesolutions tion matrix. Via Eq. (81) this again leads to Eq. (84). Alternatively, since the Green function matrix is the in- e±ika verse of the A-matrix, cf. Eq. (80), its poles are given g±(E)= , (79) by the roots of det(A) = 0. This equation is equiva- β lent to Eq. (84), as is easily shown by setting λ(+) = ±exp(−κa) in the A-matrix. for both leads. The Green function matrix in the scat- tering regioncanthen be found fromEqs.(21)and (22), which can be combined in the 3×3 matrix equation B. Fe|vacuum|Fe tunnel junction AG(E)=I, (80) As anexample of a more complex system, we consider where G (E), i,j =0,...,2 are the matrix elements of an Fe|vacuum|Fe tunnel junction where the electronic i,j G(E)andAisgivenby Eq.(74). InvertingAyieldsthe structure is treated using the local density approxima- matrix element tion of DFT.42 The calculations are based upon a tight- bindingmuffintinorbital(TB-MTO)atomicspheresap- f e2ika proximation(ASA)implementation8,19,20,21,25 ofthefor- G (E)= , (81) 2,0 2βd+(1−b)cos(ka)−ibsin(ka) malism described in Sect. III. The first step in the calculation is the self-consistent with the parameters b,d and f defined by Eq.(76). The determination of the electronic structure of the tun- relevantGreenfunctionmatrixelementfortheideallead nel junction using the layer Green function approach of is found from Eq. (34) Ref. 35. The Fe leads are oriented in the (001) direc- tion and the atoms at the Fe(001) surfaces are kept at i G(0)(E)= . (82) their unrelaxed bulk positions. For the bcc structure 0,0 2βsin(ka) and TB-MTOs,43 a principal layerin the (001) direction contains two monolayers of Fe with a thickness of 2.866 Using these results in Eq. (30) one observes that the ex- ˚A. The vacuum region is modeled by a number of such pression for the (one channel) transmission matrix ele- slices, of the same thickness, filled with “empty” spheres ment becomes identical to Eq. (75). of the same size and packing as the Fe atomic spheres. Finally one can calculate the transmission probability The atomic sphere potentials of the vacuum region and fromtheCaroliexpressiongiveninSec.VD,cf. Eq.(53). four monolayers (two principal layers) of Fe on either Using Eqs. (26), (54) and (79) one obtains side of the vacuum are calculated self-consistently while the potentials of more distant layers are kept at their ΓL =ΓR =−2βsin(ka), (83) bulk values. These potentials then form the input to a transmission calculation based on mode matching.8,20,21 for left and right leads. Using Eqs. (81), (83) and Further technical details can be found in Ref. 25. Ga0,2 = (G2,0)∗ in Eq. (53) then yields an expression for A useful quantity to extract from the self-consistent the transmissionprobabilitythatisidenticaltoEq.(77). layer calculation is the layer density of states (LDOS) It illustrates the equivalence of the different approaches ρ (E). ItisrelatedtotheretardedGreenfunctionmatrix i for calculating the transmission in this simple model. defined in Eq. (21) by In addition to providing a channel for propagating states, an impurity can also give rise to localized states, ρ (E+iη)=−π−1ImTr[G (E+iη)] (85) i i,i whose energyis outside the energyband ofthe chain, cf. Eq. (72). Such a state does not contribute to the phys- wherethetracereferstotheusuallmangularmomentum ical transmission, but the transmission amplitude has a indices characterizing MTOs. For reasons of numerical poleatanenergythatcorrespondstoalocalizedstate.41 stabilitytheretardedGreenfunctionmatrixiscalculated Within the mode matchingapproachthis correspondsto by adding a finite imaginary partto the energy; we have anenergyatwhichc becomesinfinite. Forthepresent used η = 0.0025 Ry. In the mode matching approach, S+1 modelthe energiesoflocalizedstates canbe obtainedby theLDOScanbe directlyexpressedintermsofthe wave 10 FIG.3: (Color online) Layerdensityofstates (LDOS)at the Fe(001) surface layer. Top panel: majority spin. Top curve (red, dashed): as calculated using the layer Green function method.35 Bottomcurve(blue,solid): ascalculatedusingthe mode matching approach. Both calculations use η = 0.0025 Ry. For clarity, the top curve has been displaced by 10 units along the y-axis. Bottom panel: minority spin, the LDOS is shown with a negative sign. functionsor,alternatively,intermsoftheGreenfunction matrix of Eq. (18) ρ (E)=−π−1ImTr[G′ (E)]. (86) i i,i This Green function matrix can be calculated for a real energy, but in order to make a comparisonto the results obtained with Eq. (85), we add an imaginary part, η = 0.0025 Ry. In Fig. 3 we compare the LDOS obtained using (85) and (86) for the topmost Fe monolayer of (001) Fe|vacuum|Fe where the vacuum layerwasso thick (four principal layers, corresponding to a thickness of 11.466 ˚A) that the LDOS corresponds closely to that of a free Fe(001) surface. The two curves, displaced vertically for clarity, are indistinguishable, illustrating the equiv- alence of the two Green functions defined in (18) and (21). Moreover, the LDOS are in essential agreement with results found previously for a Fe(001) surface.44 If one resolves the LDOS into contributions from dif- ferent parts of the surface Brillouin zone, then the con- tribution at Γ (k = 0) exhibits sharp peaks near the k Fermi level. These are associated with characteristic FIG. 4: (Color online) The top panels represent the band surface states found on (001) surfaces of bcc transition structures along Γ-H of bulk Fe for majority and minority metals.45 These surface states are mainly derived from spins. Panels (a-h) give the ∆1 contribution to the LDOS d3z2−r2 orbitals on the surface atoms projecting into the ρi(E) at Γ of the Fe|vacuum interface as a function of the vacuum, and they have ∆ symmetry.46,47 They become layerposition. (a)LDOSofanFelayerataninfinitedistance 1 from the interface; the double peak structure below -4 eV most clearly visible if one filters out the contribution to and single peak (of a double peak structure) above -1 eV the LDOS at Γ of the states with ∆ symmetry. 1 correspond to the ∆1 bulk bands. (b) LDOS of the 6th Fe This contribution, calculated using the mode match- layer (the surface layer is layer 1). (c) LDOS of the 3rd Fe ing approach, is shown in Fig. 4 as a function of the layer; the sharp peak in the bulk band gap corresponds to distance from the surface. In the majorityspin LDOS, a a surface state. (d) LDOS of the the surface Fe layer; the sharppeak isfoundatE −1.8eV,whichisverypromi- surface state peak has maximum amplitude. (e-h) Same for F nent in the surface Fe layer, cf. Fig. 4(d) and decays minority spin.

See more

The list of books you might like

Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.
DMCA & Copyright: Dear all, most of the website is community built, users are uploading hundred of books everyday, which makes really hard for us to identify copyrighted material, please contact us if you want any material removed.
Copyright @ PDFdrive.to, 2022 | We love our users