See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/1876527

Non Equilibrium Green's Functions for Dummies: Introduction to the One

Particle NEGF equations

Article · November 2002

Source: arXiv

CITATIONS READS
38 1,336

1 author:

Magnus Paulsson
Linnaeus University
72 PUBLICATIONS   3,425 CITATIONS   

SEE PROFILE

All content following this page was uploaded by Magnus Paulsson on 28 November 2016.

The user has requested enhancement of the downloaded file.

arXiv:cond-mat/0210519v2 [cond-mat.mes-hall] 3 Jan 2006

Non Equilibrium Green’s Functions for Dummies:

Introduction to the One Particle NEGF equations
Magnus Paulsson∗
Dept. of micro- and nano-technology, NanoDTU,
Technical University of Denmark
February 1, 2008

Abstract
Non equilibrium Green’s function methods are regularly used to calculate current and charge
densities in nanoscale (both molecular and semiconductor) conductors under bias. This method is
mainly used for ballistic conduction but may be extended to include inelastic scattering. In this tutorial
paper the NEGF equations for the current and charge density matrix are derived and explained in a
hopefully clear way.

1 Introduction
Non equilibrium Green’s function methods are regularly used to calculate current and charge densities
in nanoscale (both molecular and semiconductor) conductors under bias. An overview of the theory of
molecular electronics can be found in Ref. [1] and for semiconductor nanoscale devices see Ref. [2].
The aim of this text is to provide some intuitive explanations of one particle Green’s functions in a
compact form together with derivations of the expressions for the current and the density matrix. It is
not intended as a complete stand-alone tutorial, but rather as a complement to Ref. [1, 2, 3, 4].

2 Green’s functions
Discrete Schrödinger equation:
H|ni = E|ni (1)
We divide the Hamiltonian and wavefunction of the system into contact (H1,2 , |ψ1,2 i) and device (Hd ,
|ψd i) subspaces:     
H1 τ1 0 |ψ1 i |ψ1 i
 τ † Hd τ †   |ψd i  = E  |ψd i  (2)
1 2
0 τ2 H2 |ψ2 i |ψ2 i
∗ mpn@mic.dtu.dk

1
where τ1,2 describes the interaction between device and contacts. In general we have N contacts (H1,···,N )
connecting (τ1,···,N ) the device Hd to the reservoirs. Here we will assume that the contacts are independent,
i.e., there are no cross terms (τ ) between the different contacts.
We define the Green’s function1 :
(E − H) G(E) = I (3)

2.1 Why do we want to calculate the Green’s function?

The Green’s function gives the response of a system to a constant perturbation |vi in the Schrödinger
equation:
H|ψi = E|ψi + |vi (4)
The response to this perturbation is:

(E − H) |ψi = −|vi → (5)

|ψi = −G(E)|vi (6)

Why do we need the response to this type of perturbation? Well, it turns out that it’s usually easier
(see next section) to calculate the Green’s function than solve the whole eigenvalue problem2 and most
(all for the one-particle system) properties of the system can be calculated from the Green’s function.
E.g., the wavefunction of the contact (|ψ2 i) can be calculated if we know the wavefunction on the device
(|ψd i). From third row of Eq. 2:

H2 |ψ2 i + τ2 |ψd i = E|ψ2 i → (7)

(E − H2 ) |ψ2 i = τ2 |ψd i → (8)
|ψ2 i = g2 (E)τ2 |ψd i (9)

where g2 is the Green’s function of the isolated contact 2 ((E − H2 )g2 = I).
It is important to note that since we have an infinite system, we obtain two types of solutions for the
Green’s functions3 , the retarded and the advanced4 solutions corresponding to outgoing and incoming
waves in the contacts.
Notation: We will denote the retarded Green’s function with G and the advanced with G† (and maybe
G and GA occasionally). Here, CAPITAL G denotes the full Green’s function and its sub-matrices G1 ,
R

Gd , G1d etc. Lowercase is used for the Green’s functions of the isolated subsystems, e.g., (E − H2 ) g2 = I.
Note that by using the retarded Green’s function of the isolated contact (g2 ) in Eq. 9 we obtain the
solution corresponding to a outgoing wave in the contact. Using the advanced Green’s function (g2† ) would
give the solution corresponding to an incoming wave.
1 Others may (and do) use the opposite sign.
2 Especially for infinite systems.
3 When the energy coincides with energy band of the contacts there are two solutions corresponding to outgoing or

incoming waves in the contacts.

4 In practice these two solutions are usually obtained by adding an imaginary part to the energy. By taking the limit to

zero of the imaginary part one of the two solutions is obtained. If the limit → 0+ is taken the retarded solution is found,
→ 0− gives the advanced. This can be seen from the Fourier transform of the time dependent Green’s function.

2
2.2 Self-Energy
The reason for calculating the Green’s function is that it is easier that solving the Schrödinger equation.
Also, the Green’s function of the device (Gd ) can be calculated separately without calculating the whole
Green’s function (G). From the definition of the Green’s function we obtain:
    
E − H1 −τ1 0 G1 G1d G12 I 0 0
 −τ † E − Hd −τ2†   Gd1 Gd Gd2  =  0 I 0  (10)
1
0 −τ2 E − H2 G21 G2d G2 0 0 I

Selecting the three equations in the second column:

(E − H1 ) G1d − τ1 Gd = 0 (11)
−τ1† G1d + (E − Hd ) Gd − τ2† G2d = I (12)
(E − H2 ) G2d − τ2 Gd = 0 (13)

We can solve Eqs. 11 and 13 for G1d and G2d :

G1d = g1 τ1 Gd (14)
G2d = g2 τ2 Gd (15)

substitution into Eq. 12 gives:

− τ1† g1 τ1 Gd + (E − Hd ) Gd − τ2† g2 τ2 Gd = I (16)

from which Gd is simple to find:

−1
Gd = (E − Hd − Σ1 − Σ2 ) (17)
where Σ1 = τ1† g1 τ1and Σ2 = τ2† g2 τ2are the so called self-energies.
Loosely one can say that the effect of the contacts on the device is to add the self-energies to the
device Hamiltonian since when we calculate the Green’s function on the device we just calculate the
Green’s function for the effective Hamiltonian Heffective = Hd + Σ1 + Σ2 . However, we should keep in mind
that we can only do this when we calculate the Green’s function. The eigen-values and -vectors of this
effective Hamiltonian are not quantities we can interpret easily.
For “normal” contacts, the surface Green’s functions g1 and g2 used to calculate the self-energies are
usually calculated using the periodicity of the contacts, this method is described in detail in appendix B
of Ref. [3] and in section 3 of Ref. [2].

2.3 The spectral function

Another important use of the Green’s function is the spectral function:

A = i G − G†


(18)

which gives the DOS and all solutions to the Schrödinger equation!
To see this we first note that for any perturbation |vi we get two solutions (|ψ R i and |ψ A i) to the
perturbed Schrödinger equation:
(E − H)|ψi = −|vi (19)

3
from the advanced and retarded Green’s functions:

|ψiR = −G|vi (20)

|ψiA = −G† |vi (21)

The difference of these solutions (|ψ R i − |ψ A i = (G − G† )|vi = −iA|vi) is a solution to the Schrödinger
equation:
(E − H)(|ψ R i − |ψ A i) = (E − H)(G − G† )|vi = (I − I)|vi = 0 (22)
which means that |ψi = A|vi is a solution to the Schrödinger equation for any vector |vi!
To show that the spectral function actually gives all solutions to the Schrödinger equation is a little
bit more complicated and we need the expansion of the Green’s function in the eigenbasis:

1 X |kihk|
G= = (23)
E + iδ − H E + iδ − ǫk
k

where the δ is the small imaginary part (see footnote 3), |ki’s are all eigenvectors5 to H with the corre-
sponding eigenvalues ǫk . Expanding the spectral function in the eigenbasis gives:
 
1 1
A = i − (24)
E + iδ − H E − iδ − H
 
X 1 1
= i |kihk| − (25)
E + iδ − ǫk E − iδ − ǫk
k
X 2δ
= |kihk| 2 (26)
k
(E − ǫk ) + δ 2

where δ is our infinitesimal imaginary part of the energy. Letting δ go to zero gives:
X
A = 2π δ (E − ǫk ) |kihk| (27)
k

(here δ(E − ǫk ) is the delta function) which can be seen since (E−ǫ2δ)2 +δ2 goes to zero everywhere but at
k
E = ǫk , integrating over E (with a test function) gives the 2πδ (E − ǫk ) factor. Eq. 27 shows that the
spectral function gives us all solutions to the Schrödinger equation.

3 Response to an incoming wave

In the non-equilibrium case, reservoirs with different chemical potentials will inject electrons and occupy
the states corresponding to incoming waves in the contacts. Therefore, we want to find the solutions
corresponding to these incoming waves.
Consider contact 1 isolated from the other contacts and the device. At a certain energy we have
solutions corresponding to an incoming wave that is totally reflected at the end of the contact. We will
denote these solutions with |ψ1,n i where 1 is the contact number and n is a quantum number (we may
5 Normalized!

4
have several modes in the contacts). We can find all these solutions from the spectral function a1 of the
isolated contact (as described above).
Connecting the contacts to the device we can calculate the wavefunction on the whole system caused by
the incoming wave in contact 1. To do this we note that a wavefunction should be of the form |ψ1,n i+|ψ R i
where |ψ1,n i is the totally reflected wave and |ψ R i is the retarded response of the whole system. Putting
in the anzats |ψ1,n i + |ψ R i into the Schrödinger equation gives:
 
H1 + τ1 +
 Hd + τ † + τ † +  |ψ1,n i + |ψ R i = E |ψ1,n i + |ψ R i



1 2 (28)
H2 + τ2
   
E|ψ1,n i+ H1 + τ1 +
 τ † |ψ1,n i+  +  Hd + τ † + τ † +  |ψ R i = E |ψ1,n i + |ψ R i


1 1 2 (29)
0 H2 + τ2
H|ψ R i = E|ψ R i − τ1† |ψ1,n i (30)

(note the slight change in notation) and we see that |ψ R i is noting else the response of the whole system
to a perturbation of −τ1† |ψ1,n i, c.f., Eq. 9:

|ψ R i = Gτ1† |ψ1,n i (31)

It is important to realize that the scattering states generated from Eq. 31, using all possible incoming
waves from each contact, form a complete6 ON set of solutions to the full Schrödinger equation[5]. Note
that we have chosen the retarded response which means that the only part of the wave that is traveling
towards the device is the incoming wave (part of |ψ1,n i). We will make full use of this fact below.
It will be useful to have the expressions for the device wavefunction |ψd i and contact wavefunction
(|ψ1,2 i). The device part is straightforward:

|ψd i = Gd τ1† |ψ1,n i (32)

and from Eq. 9 or Eq. 15:

|ψ2 i = g2 τ2 |ψd i = g2 τ2 Gd τ1† |ψ1,n i (33)
Note, to calculate the wavefunction in the contact containing the incoming wave (contact 1) we need to
add the incoming wave, giving a slightly more complicated expression:
 
|ψ1 i = 1 + g1 τ1 Gd τ1† |ψ1,n i (34)

Knowing the wavefunctions corresponding to incoming waves in different contacts enables us to fill up
the different solutions according to the electron reservoirs filling the contacts.

4 Charge density matrix

In the non equilibrium case we are often interested in two quantities: the current and the charge density
matrix. Lets start with the charge density (which allows us to use a self-consistent scheme to describe
charging).
6 Except for localized states in the device region.

5
 The charge density matrix is defined as:
X
ρ= f (k, µ)|ψk ihψk | (35)
k

where the sum runs over all states with the occupation number f (Ek , µ) (pure density matrix) (note the
similarity with the spectral function A, in equilibrium you find the density matrix from A and not as
described below). In our case, the occupation number is determined by the reservoirs filling the incoming
waves in the contacts such that:
1
f (Ek , µ1 ) = (36)
1 + e −µ1 )/kB T
(E k

is the Fermi-Dirac function with the chemical potential (µ1 ) and temperature (T ) of the reservoir respon-
sible for injecting the electrons into the contacts.
The wavefunction on the device given by an incoming wave in contact 1 (see Eq. 32) is:

|ψd,k i = Gd τ1† |ψ1,k i (37)

Adding up all states from contact 1 gives:

Z∞ X
ρd [contact 1] = dE f (E, µ1 )δ(E − Ek )|ψd,k ihψd,k | (38)
E=−∞ k

Z∞ X
= dE f (E, µ1 ) δ(E − Ek )Gd τ1† |ψ1,k ihψ1,k |τ1 G†d (39)
E=−∞ k

Z∞ X
!
= dE f (E, µ1 )Gd τ1† δ(E − Ek )|ψ1,k ihψ1,k | τ1 G†d (40)
E=−∞ k

Z∞
a1
= [Eq. 27] = dE f (E, µ1 )Gd τ1† τ1 G†d (41)
2π
E=−∞
 
introducing the new quantity Γ1 = τ1† a1 τ1 = i Σ1 − Σ†1 we obtain the simple formula:

Z∞
1
ρ[from contact 1] = dE f (E, µ1 )Gd Γ1 G†d (42)
2π
E=−∞

The total charge density thus becomes a sum over all contacts:
Z∞
2 (for spin) X
ρ= dE f (E, µi )Gd Γi G†d (43)
2π i
E=−∞

6
5 Probability Current
Having different chemical potentials in the reservoirs filling the contacts gives rise to a current. In the
next section we will calculate this current in a similar way as the charge density was calculated. But to
do this we need an expression for the current from the wavefunction.
In the continuum case we can calculate the current from the velocity operator. However, for a discrete
Hamiltonian it is not so clear what the velocity operator is. Therefore, we derive an expression for the
current from the continuity equation (using two contacts). In steady-state, the probability to find an
electron on the device ( |ψi |2 where the sum runs over the device subspace) is conserved:
P
i

|ψi |2
P
∂ X ∂hψ|iihi|ψi X  ∂hψ|ii 
i ∂hi|ψi
0 = = = hi|ψi + hψ|ii (44)
∂t i
∂t i
∂t ∂t
i X i
= (hψ|H|iihi|ψi − hψ|iihi|H|ψi) = (hψ|H|ψd i − hψd |H|ψi) (45)
h̄ i h̄
i  
= hψ|Hd + τ1 + τ2 |ψd i − hψd |Hd + τ1† + τ2† |ψi (46)
h̄
i h i h i
= hψ1 |τ1 |ψd i − hψd |τ1† |ψ1 i + hψ2 |τ2 |ψd i − hψd |τ2† |ψ2 i (47)
h̄
We interpret the term in the first (square) bracket as the incoming probability current into the device
from contact 1 and the second bracket from contact 2. Generalizing to an arbitrary contact j gives us the
electric current (at one energy) as the charge (-e) times the probability current:
ie  
ij = − hψj |τj |ψd i − hψd |τj† |ψj i (48)
h̄
where ij is defined as positive for a current from the contacts into the device. We can now put in the
expressions for the wavefunctions in the same way as for the density matrix.

6 Electrical Current
To calculate the total current through the device we only need to put in the wavefunction of the device
and the contacts (|ψd i, |ψ1 i, |ψ2 i) from Eqs. 32, 34 and 9 and add all the contributions together. Thus the
current into the device from a incoming wave of one energy (E) in contact 1 (|ψ1,n i) through the coupling
defined by τ2 is:
ie  
i2 from 1 = − hψ2 |τ2 |ψd i − hψd |τ2† |ψ2 i (49)
h̄
ie
= − (hψ1,n |τ1 G†d τ2† g2† τ2 Gd τ1† |ψ1,n i − hψ1,n |τ1 G†d τ2† g2 τ2 Gd τ1† |ψ1,n i) (50)
h̄
ie  
= − hψ1,n |τ1 G†d τ2† g2† − g2 τ2 Gd τ1† |ψ1,n i) (51)
h̄
e
= hψ1,n |τ1 G†d Γ2 Gd τ1† |ψ1,n i (52)
h̄

7
Adding over the modes n and noting that the levels are filled from the reservoir connected to contact 1
gives (2 for spin):
Z∞
e X
I2 from 1 = 2 dE f (E, µ1 ) δ(E − En )hψ1,n |τ1 G†d Γ2 Gd τ1† |ψ1,n i (53)
h̄ n
E=−∞
Z∞
2e X
= dE f (E, µ1 ) δ(E − En )hψ1,n |τ1 |mihm|G†d Γ2 Gd τ1† |ψ1,n i (54)
h̄ m,n
E=−∞
Z∞ !
2e X X
= f (E, µ1 ) hm|G†d Γ2 Gd τ1† δ(E − En )|ψ1,n ihψ1,n | τ1 |mi (55)
h̄ m n
E=−∞
Z∞
2e X a1
= dE f (E, µ1 ) hm|G†d Γ2 Gd τ1† τ1 |mi (56)
h̄ m
2π
E=−∞
Z∞
e  
= dE f (E, µ1 )Tr G†d Γ2 Gd Γ1 (57)
πh̄
E=−∞

To get the total current through the device the current from contact two have to be subtracted away:
Z∞
e  
I= dE (f (E, µ1 ) − f (E, µ2 )) Tr G†d Γ2 Gd Γ1 (58)
πh̄
E=−∞

which is exactly the Landauer formula for the current.

Acknowledgments
The author is grateful for the inspiration from Supriyo Datta and Mads Brandbyge. The revised version
was initiated from discussions with Casper Krag and Carsten Rostgård.

References
[1] M. Paulsson, F. Zahid, and S. Datta. Nanoscience, Engineering and Technology Handbook, chapter
Resistance of a Molecule. CRC Press, 2002. Editors, D. Brenner, S. Lyshevski and G. Iafrate, cond-
mat/0208183.
[2] P. S. Damle, A. W. Ghosh, and S. Datta. Molecular nanoelectronics, chapter Theory of nanoscale
device modeling. 2002. Editor M. Reed.
[3] F. Zahid, M. Paulsson, and S. Datta. Advanced Semiconductors and Organic Nano-Techniques, chapter
Electrical Conduction through Molecules. Academic press, 2002. Editor H. Markoc.
[4] S. Datta. Electronic transport in mesoscopic systems. Cambridge University Press, Cambridge, UK,
1997.

8
 [5] L. E. Ballentine. Quantum mechanics. World Scientific, 1998. see chapter 16, ”Scattering”.

View publication stats