Seminario J.M. Molecular and Nano Electronics. Analysis, Design and Simulation
Подождите немного. Документ загружается.
254 G. Stefanucci et al.
Due to the function, we have t>t
for c
R
= 0. In the second term on the r.h.s.
b
A
¯
tt
contains a t
−
¯
t and hence it must be t>
¯
t; therefore we can replace the
difference in the square bracket with a
R
. Then we break the first term on the r.h.s. into
two pieces by inserting functions: one for
¯
t<t
and the other for
¯
t>t
. In compact
notation, we end up with
c
R
=a
R
·b
R
c
A
=a
A
·b
A
(17)
where the second relation can be proven in a similar way. It is worth noting that in the
expressions for c
R
and c
A
no integration along the imaginary track is required. For later
purposes we also define the Matsubara function k
M
with both the arguments in
the interval 0 −i:
k
M
≡ kz = z
=
(18)
It is straightforward to prove the following identities
c
=a
R
·b
+a
∗b
M
c
=a
·b
A
+a
M
∗b
c
M
=a
M
∗b
M
(19)
Finally, we consider the case of a Keldysh function kz z
multiplied on the left
by a scalar function lz. The function k
l
z z
= lzkz z
= z z
lzk
>
z z
+
z
zlzk
<
z z
and hence belongs to the Keldysh space. The Keldysh components
can be extracted using the definitions (10, 11, 14, 18). Choosing for instance z = t
+
and z
=t
−
we find k
>
l
t t
=ltk
>
t t
and similarly for z =t
−
and z
=t
+
we find
k
<
l
t t
= ltk
<
t t
. More generally, the function l is simply a prefactor: k
x
l
= lk
x
,
where x is one of the Keldysh components (≶,R A, , M). The same is true for
k
r
z z
= kz z
rz
, where rz
is a scalar function: k
x
r
=k
x
r.
3. Quantum transport using TDDFT and NEGF
3.1. Merging the Keldysh and TDDFT formalisms
The one-particle scheme of TDDFT corresponds to a fictitious system of non-
interacting electrons described by the Kohn–Sham (KS) Hamiltonian
ˆ
H
s
z =
drdr
†
rrH
s
zr
r
. The potential v
s
rt experienced by the electrons in
the free-electron Hamiltonian H
s
t is called the KS potential and it is given by the sum
of the external potential, the Coulomb potential of the nuclei, the Hartree potential and
the exchange-correlation potential v
xc
. The latter accounts for the complicated many-
body effects and is obtained from an exchange-correlation action functional, v
xc
rt=
A
xc
n/nrt (as pointed out in [17], the causality and symmetry properties require
that the action functional A
xc
n is defined on the Keldysh contour). A
xc
is a functional
of the density and of the initial density matrix. In our case, the initial density matrix is
the thermal density matrix which, due to the extension of the Hohenberg–Kohn theorem
[18] to finite temperatures [19], also is a functional of the density. We should mention
here that an alternative formulation based on TDDFT has been recently proposed by
Di Ventra and Todorov [20]. In their approach the system is initially unbalanced and
Time-dependent transport phenomena 255
therefore the exchange-correlation functional depends on the initial state and not only
on the density.
The fictitious Keldysh–Green function z z
of the KS system satisfies a one-
particle equation of motion
i
−→
d
dz
1−H
s
z
z z
= 1z −z
z z
−i
←−
d
dz
1−H
s
z
=1z −z
(20)
with KMS boundary conditions (9). In Eqs (20) the arrow specifies where the derivative
along the contour acts. The left and right equations of motion are equations on the
extended Keldysh contour of Figure 2c and z −z
=
d
dz
z z
=−
d
dz
z z
. For
any z = z
, the equations of motion are solved by the evolution operator on the con-
tour Sz z
= T
K
e
−i
z
z
d¯zH
s
¯z
, since i
−→
d
dz
Sz z
= H
s
zSz z
and Sz z
−i
←−
d
dz
=
Sz z
H
s
z
. Therefore, any Green function
z z
= z z
Sz 0
−
f
>
S0
−
z
+z
zSz 0
−
f
<
S0
−
z
(21)
with f
≶
constrained by f
>
−f
<
=−i1, is a solution of Eqs (20). In order to fix the
matrix f
>
or f
<
we impose the KMS boundary conditions. The matrix H
s
z equals
H
s
for any z on the vertical track, meaning that S−i 0
−
= e
−H
s
. Equations (9)
then implies f
<
=−e
−
H
s
−
f
>
, and taking into account the constraint f
>
−f
<
=−i1
we conclude that f
<
= ifH
s
, where f = 1/e
−
+1 is the Fermi distribution
function. The matrix f
>
takes the form f
>
=ifH
s
−1.
The Green function z z
for a system of noninteracting electrons is now completely
fixed. Let us consider some Keldysh–Green functions. For z =t
+
and z
=t
−
we have
the greater Green function while for z = t
−
and z
= t
+
we have the lesser Green
function
>
t t
= iSt 0fH
s
−1S0t
<
t t
= iSt 0fH
s
S0t
(22)
Both
>
and
<
depend on the initial distribution function fH
s
. The diagonal
matrix element of −i
<
is nothing but the time-dependent value of the local electron
density nrt, while i
>
gives the local density of holes. Another way of writing
−i
<
is in terms of the eigenstates
s
0 of H
s
with eigenvalues
s
. From the
time-evolved eigenstate
s
t=St 0
s
0, we can calculate the time-dependent
wavefunction
s
rt=r
s
t. Inserting
s
s
0
s
0 in the expression for
<
we find −ir
<
t t
r
=
s
f
s
s
r t
∗
s
r
t
, which for t = t
reduces to the
time-dependent density matrix. Knowing the greater and lesser Green functions we can
also calculate
RA
. Taking into account the definitions (14) we find
R
t t
=−it −t
St t
A
t t
= it
−tSt t
=
R
t
t
†
(23)
In the above expressions for
RA
, the Fermi distribution function has disappeared.
The information carried by
RA
is the same as contained in the one-particle evolution
256 G. Stefanucci et al.
operator. There is no information on how the system is prepared (how many particles,
how they are distributed, etc.). We use this observation to rewrite
≶
in terms of
RA
≶
t t
=
R
t 0
≶
0 0
R
0t
(24)
Thus,
≶
is completely known once we know how to propagate the one-electron
orbitals in time and how they are populated before the system is perturbed [4, 21]. For
later purposes, we also observe that an analogous relation holds for
t = i
R
t 0
0
t =−i
0
A
0 t (25)
3.2. Total current using TDDFT
The fictitious of the KS system will not, in general, give correct one-particle properties.
However, by definition,
<
gives the correct density nrt=−ir
<
t tr. Also
total currents are correctly given by TDDFT. If, for instance, I
is the total current from
a particular region we have
I
t =−e
dr
d
dt
nrt= e
dr i
d
dt
r
<
t tr (26)
where the space integral extends over the region (e is the electron charge). We stress
here that I
is the electronic current (the direction of the current coincides with the
direction of the electrons).
At this point, it is convenient to partition the system into three main regions: a
central region C consisting of the junction and a few atomic layers of the left and right
electrodes and two regions L and R which describe the left and right bulk electrodes.
According to this partitioning, the KS Hamiltonian H
s
can be written as a 3 ×3 block
matrix, and the left equation of motion in (20) reads
⎧
⎨
⎩
i
d
dz
1−
⎡
⎣
H
LL
z H
LC
0
H
CL
H
CC
z H
CR
0 H
RC
H
RR
z
⎤
⎦
⎫
⎬
⎭
z z
= z −z
1 (27)
with
z z
=
⎡
⎣
LL
z z
LC
z z
LR
z z
CL
z z
CC
z z
CR
z z
RL
z z
RC
z z
RR
z z
⎤
⎦
(28)
(a similar equation is easily obtained for the right equation of motion). Choosing z
on the forward branch of the Keldysh contour and z
on the backward branch of the
same contour, we obtain a left and a right equation for the lesser Green function. These
equations can be used to get rid of the time derivative in Eq. (26). We find for =L R
I
t = e
dr ri
d
dt
<
t tr
=e
dr rH
C
<
C
t t −
<
C
t tH
C
r=2e Re
Tr
C
Q
t
(29)
Time-dependent transport phenomena 257
where
Q
t ≡
<
C
t tH
C
=
R
t 0
<
0 0
A
0t
C
H
C
=
R
CC
t 0
<
CC
0 0
A
C
0tH
C
+
=LR
R
C
t 0
<
C
0 0
A
C
0tH
C
+
=LR
R
CC
t 0
<
C
0 0
A
0tH
C
+
=LR
R
C
t 0
<
0 0
A
0tH
C
(30)
is a one-particle operator in the central region C and Tr
C
denotes the trace over a
complete set of one-particle states of C. Let us express the quantity Q
in terms of
the Green function
CC
projected in the central region. We introduce the uncontacted
Green function g which obeys Eqs (20) with H
C
=H
C
=0,
⎧
⎨
⎩
i
d
dz
1−
⎡
⎣
H
LL
z 00
0 H
CC
z 0
00H
RR
z
⎤
⎦
⎫
⎬
⎭
gz z
= z −z
1 (31)
where
gz z
=
⎡
⎣
g
LL
z z
00
0 g
CC
z z
0
00g
RR
z z
⎤
⎦
(32)
and the same KMS boundary conditions as . The uncontacted g allows us to convert
Eqs (20) into an integral equation which entails the KMS boundary conditions for :
z z
= gz z
+
d¯z gz ¯zH
off
¯z z
=gz z
+
d¯z z ¯zH
off
g¯z z
(33)
being the extended Keldysh contour of Figure 2c and H
off
is the off-diagonal part of
H
s
. Using the relations (17) of Section 2.3 we find
RA
C
=
RA
CC
·H
C
g
RA
A
=
g
A
+g
A
H
C
·
A
CC
·H
C
g
A
(34)
In Eq. (30) all matrix elements of
<
are evaluated at times 0 0. From Eq. (15) we
see that c
<
0 0 =
a
∗b
0 0, due to the functions in the retarded and advanced
components. Therefore
<
C
0 0 =
g
H
C
∗
CC
0 0
<
C
0 0 =
CC
∗H
C
g
0 0 (35)
and exploiting the first two relations in Eq. (19) we also find that
<
0 0 =
g
<
0 0 +
g
H
C
∗
M
CC
∗H
C
g
0 0 (36)
258 G. Stefanucci et al.
Substituting Eqs (34–36) into Eq. (30) and using the identities (24 and 25) for the
Green function g, we obtain the following expression for Q
t:
Q
t =
=LR
G
R
·
<
·
+G
A
·
A
t t
+
=LR
G
R
·
∗G
M
∗
·
+G
A
·
A
t t
+i
=LR
G
R
t 0
G
∗
·
+G
A
·
A
0t
+
G
R
t 0G
<
0 0 −i
G
R
·
∗G
t 0
G
A
·
A
0 t (37)
where we have used the short-hand notation G ≡
CC
and
z z
=
=LR
z z
= H
C
g
z z
H
C
(38)
is the so-called embedding self-energy which accounts for hopping in and out of
region C.
Having the quantity Q
t we can calculate the exact total current I
t of an interact-
ing system of electrons. Eq. (29) allows for studying transient effects and more generally
any kind of time-dependent current responses. In the long time limit
lim
t→
Q
t =
G
R
·
<
+G
R
·
<
·G
A
·
A
t t (39)
provided G and tend to zero when the separation between their time arguments
increases (in this case, it is only the first term on the r.h.s. of Eq. (37) that does not
vanish). This condition is not stringent and is fulfilled provided the electrode states
form a continuum and the local density of states in the central region C is a smooth
function. In the next section we investigate under what circumstances a steady current
I
develops in the long-time limit. We will also discuss the dependence of I
on the
history of the external potential.
3.3. Steady state and history dependence
In this section we show that a steady state develops provided (1) the KS Hamiltonian
H
s
t globally converges to an asymptotic KS Hamiltonian H
s
when t →and (2)
the electrodes form a continuum of states (thermodynamic limit), and the local density
of states is a smooth function in the central region.
Let us define the asymptotic KS Hamiltonian of electrode as H
=lim
t→
H
t.
The retarded/advanced component of the uncontacted Green function g behaves like
lim
t→
g
R
t 0 = ie
−iH
t
lim
t→
g
A
0t=−i
†
e
iH
t
(40)
where is a unitary operator and it is defined according to
= lim
t→
T
e
−i
t
0
H
t
dt
e
−iH
t
(41)
Time-dependent transport phenomena 259
T being the time-ordering operator. In terms of diagonalizing one-body states
m
of
H
with eigenvalues
m
, the lesser component of the embedding self-energy, defined
in Eq. (38), can be written as
lim
tt
→
<
t t
= lim
tt
→
H
C
g
R
t 0g
<
0 0g
A
0t
H
C
=i
mm
e
−i
m
t−
m
t
×H
C
m
m
fH
0
†
m
m
H
C
(42)
where we have taken into account that g
<
0 0 = ifH
0. The left and right
contraction with a nonsingular hopping matrix H
C
causes a perfect destructive inter-
ference for states with
m
−
m
1/t +t
and hence the restoration of translational
invariance in time
lim
tt
→
<
t t
= i
m
f
m
m
e
−i
m
t−t
(43)
where f
m
=
m
fH
0
†
m
while
m
=H
C
m
m
H
C
. In principle,
there may be degeneracies which require a diagonalization to be performed for states on
the energy shell. The above dephasing mechanism is the key ingredient for a steady state
to develop. Substituting Eq. (43) into Eq. (39), we obtain for the steady state current
I
S
=−2e
m
f
m
Tr
C
G
R
m
m
G
A
m
Im
A
m
−2e
m
f
m
Tr
C
m
ImG
R
m
(44)
with
G
RA
=
1
1
C
−H
CC
−
RA
(45)
The imaginary part of G
R
is simply given by G
R
Im
R
G
A
. By definition, we have
RA
= H
C
1
1
−H
±i
H
C
(46)
and hence
Im
RA
=∓
m
−
m
m
(47)
Using the above identity, the steady-state current can be rewritten in a Landauer-like
[22] form
I
S
R
=−e
m
f
mL
mL
−f
mR
mR
=−I
S
L
(48)
260 G. Stefanucci et al.
In the above formula
mR
=
n
nL
mR
and
mL
=
n
nR
mL
are the TDDFT transmission
coefficients expressed in terms of the quantities
n
m
=2
m
−
n
Tr
C
G
R
m
m
G
A
n
n
=
m
n
(49)
Despite the formal analogy with the Landauer formula, Eq. (48) contains an important
conceptual difference since f
m
is not simply given by the Fermi distribution function.
For example, if the induced change in effective potential varies widely in space deep
inside the electrodes, the band structure of the -electrode Hamiltonian H
0
†
might differ from that of H
. However, for metallic electrodes with a macroscopic
cross section the switching on of an electric field excites plasmon oscillations, which
dynamically screen the external disturbance. Such a metallic screening prevents any
rearrangements of the initial equilibrium bulk-density, provided the time-dependent
perturbation is slowly varying during a typical plasmon time-scale (which is usually less
than a femtosecond). Thus, the KS potential v
s
undergoes a uniform time-dependent shift
deep inside the left and right electrodes and the KS potential-drop is entirely limited to
the central region. Denoting with v
t the difference in electrode between the KS
potential at time t and the KS potential at negative times, v
t =v
s
r ∈ t −v
s
r ∈
0, to leading order in 1/N we then have
H
t = H
0 +1
v
t (50)
meaning that H
=H
0 +1
v
. Hence, except for corrections which are of lower
order with respect to the system size, H
0
†
=H
0 and
f
m
=f
m
−v
(51)
The formula for the current can be further manipulated when Eq. (51) holds. Let us
write the embedding self-energy as the sum of a real and imaginary part
RA
=
∓i
/2. Using Eq. (47) we can rewrite the transmission coefficients as
mR
=Tr
C
G
R
mR
mR
G
A
mR
L
mR
(52)
mL
=Tr
C
G
R
mL
mL
G
A
mL
R
mL
(53)
Substituting these expressions in Eq. (48) and taking into account Eq. (51) we obtain
I
S
R
=−e
d
2
f −v
L
−f −v
R
Tr
C
G
R
L
G
A
R
(54)
In the above equation the Green functions are calculated from Eq. (45). The
Hamiltonian H
CC
is the KS Hamiltonian H
s
t → projected on region C and can be
obtained by evolving the system for very long times. According to the Runge–Gross
theorem, H
CC
depends on how the system was prepared at t =0 (in our case the system
is contacted and in thermal equilibrium) and on the full history of the time-dependent
density. Therefore, the use of Eq. (54) in the context of static DFT is generally not
correct. Indeed, static DFT is an equilibrium theory while here we are dealing with a
Time-dependent transport phenomena 261
nonequilibrium process. One might argue that in the linear-response regime, the static
DFT approach is free from the above criticism. Unfortunately, this is not the case.
Denoting with v
the small change v
of the effective potential in electrode and
with I
S
R
the corresponding current response, to first order in v
Eq. (54) yields
I
S
R
=−e
d
2
f
Tr
C
G
R
0
0L
G
A
0
0R
v
R
−v
L
(55)
The Green functions and the ’s in Eq. (55) refer to the system in equilibrium and
static DFT approaches can be used to evaluate the trace. However, DFT is not enough
to calculate the change v
. Indeed
v
= lim
t→
lim
x→±
v
ext
rt+V
H
rt+v
xc
rt
(56)
where x is the longitudinal coordinate, the plus sign applies for = R and the minus
sign for =L. In the above equation v
ext
is the external potential and V
H
is the Hartree
potential; their sum gives the electrostatic Coulomb potential v
C
,
v
C
= lim
t→
lim
x→±
v
ext
rt+V
H
rt
(57)
The variation v
xc
of the exchange-correlation potential can be expressed in terms of
the exchange-correlation kernel f
xc
rtr
t
= v
xc
r t/nr
t
v
xc
= lim
t→
lim
x→±
v
xc
rt= lim
t→
lim
x→±
dr
dt
f
xc
rtr
t
nr
t
(58)
The kernel f
xc
depends only on the difference t −t
. We denote by f
xc
r
the
Fourier transform of f
xc
evaluated at x =±for =R L. Then
v
xc
= lim
t→
d
2
e
−it
dr
f
xc
r
nr
(59)
with nr the Fourier transform of nrt. Rewriting v
as v
C
+v
xc
and taking into account Eq. (59), the current response I
S
R
in Eq. (55) can also be
written as
I
S
R
=−e
d
2
f
T
v
RC
−v
LC
+lim
t→
d
2
e
−it
×
dr
f
Rxc
r
−f
Lxc
r
nr
(60)
with T = Tr
C
G
R
0
0L
G
A
0
0R
. At zero temperature f/ =
−
F
, with
F
being the Fermi energy, and Eq. (60) becomes
I
S
R
=G
KS
F
v
RC
−v
LC
+lim
t→
d
2
e
−it
×
dr
f
Rxc
r
−f
Lxc
r
nr
(61)
262 G. Stefanucci et al.
where G
KS
F
=−eT
F
/2 is the conductance of the KS system. We conclude
that also in the linear-response regime static DFT is not appropriate for calculating
the conductance since dynamical exchange-correlation effects might contribute through
the last term in Eq. (61). Equation (61) can also be obtained within the framework of
time-dependent current density functional theory as it has been shown in [23].
We emphasize that the steady-state current in Eq. (48) results from a pure dephasing
mechanism in the fictitious noninteracting problem. The damping effects of scattering
are described by A
xc
and v
xc
. Furthermore, the current depends only on the asymptotic
value of the KS potential, v
s
rt →. However, v
s
rt → might depend on the
history of the external applied potential and the resulting steady-state current might be
history dependent. In these cases the full-time evolution cannot be avoided. In the case
of Time Dependent Local Density Approximation (TDLDA), the exchange-correlation
potential v
xc
depends only locally on the instantaneous density and has no memory at all.
If the density tends to a constant, so does the KS potential v
s
, which again implies that
the density tends to a constant. Owing to the nonlinearity of the problem there might
still be more than one steady-state solution or none at all. We are currently investigating
the possibility of having more than one steady-state solution.
4. Quantum transport: A practical scheme based on TDDFT
The theory presented in the previous sections allows us to calculate the time-dependent
current in terms of the Green function
CC
= G projected in the central region. In
practice, it is computationally very expensive to propagate Gz z
in time (because it
depends on two time variables) and also calculate Q
from Eq. (37). Here we describe
a feasible numerical scheme based on the propagation of KS orbitals. We remind the
reader that our electrode–junction–electrode system is infinite and non-periodic. Since
one can, in practice, only deal with finite systems we will propagate KS orbitals projected
in the central region C by applying the correct boundary conditions [10].
We specialize the discussion to nonmagnetic systems at zero temperature and we
denote with
s
r 0 ≡r
s
0 the eigenstates of H
s
t < 0. The time dependent
density can be computed in the usual way by nrt=
occ
s
rt
2
, where the sum
is over the occupied Kohn–Sham orbitals and
s
t is the solution of the KS equation
of TDDFT i
d
dt
s
t=H
s
t
s
t. Using the continuity equation, we can write the
total current I
t of Eq. (26) as
I
t =−e
occ
dr ·Im
∗
s
rt
s
rt
=−e
occ
S
d
ˆ
n ·Im
∗
s
rt
s
rt
(62)
where
ˆ
n is the unit vector perpendicular to the surface element d and the surface
S
is perpendicular to the longitudinal geometry of our system. From Eq. (62) we
conclude that in order to calculate I
t we only need to know the time-evolved KS
orbitals in region C. This is possible provided we know the dynamics of the remote
parts of the system. As at the end of Section 3.3, we restrict ourselves to metallic
electrodes. Then, the external potential and the disturbance introduced by the device
Time-dependent transport phenomena 263
region are screened deep inside the electrodes. As the system size increases, the remote
parts are less disturbed by the junction and the density inside the electrodes approaches
the equilibrium bulk-density. Thus, the macroscopic size of the electrodes leads to an
enormous simplification since the initial-state self-consistency is not disturbed far away
from the constriction. Partitioning the KS Hamiltonian as in Eq. (27), the time-dependent
Schrödinger equation reads
i
d
dt
⎡
⎣
L
C
R
⎤
⎦
=
⎡
⎣
H
LL
H
LC
0
H
CL
H
CC
H
CR
0 H
RC
H
RR
⎤
⎦
⎡
⎣
L
C
R
⎤
⎦
(63)
where
is the projected wave-function onto the region =L R C. We can solve
the differential equation for
L
and
R
in terms of the retarded Green function g
R
.
Then, we have for =L R
t=ig
R
t 0
0+
t
0
dt
g
R
t t
H
C
C
t
(64)
Using Eq. (64), the equation for
C
can be written as
i
d
dt
C
t=H
CC
t
C
t+
t
0
dt
R
t t
C
t
+i
=LR
H
C
g
R
t 0
0 (65)
where
R
=
=LR
H
C
g
R
H
C
, in accordance with Eq. (38). Thus, for any given
KS orbital we can evolve its projection onto the central region by solving Eq. (65) in
region C. Equation (65) has also been derived elsewhere (for static Hamiltonians) [24].
To summarize, all the complexity of the infinite electrode–junction–electrode system
has been reduced to the solution of an open quantum-mechanical system (the central
region C) with proper time-dependent boundary conditions.
Equation (65) is the central equation of our numerical approach to time-dependent
transport. It is a reformulation of the original time-dependent Schrödinger equation (63)
of the full system in terms of an equation for the central (device) region only. The
coupling to the leads is taken into account by the lead Green functions g
R
= L R.
Equation (65) has the structure of a time-dependent Schrödinger equation with two extra
terms. The first term describes the injection of particles induced by a nonvanishing
projection of the initial wave-function onto the leads. The second term involves the
self-energy
R
and the wavefunction in the central region at previous times during the
propagation. We will denote it as the memory integral. We should remark here that
these memory effects are of different origin than those which are usually discussed
in the context of TDDFT [25, 26]. The latter ones arise from the dependence of the
exchange-correlation functional on the full history of the time-dependent density. Most
density-based functionals used at present rely on the adiabatic approximation therefore
ignoring the functional dependence on past time-dependent densities [27].
Equation (65) is first order in time, therefore we need to specify an initial state which
is to be propagated. We want to study the time evolution of systems perturbed out of
their equilibrium ground state. Of course, the ground state of our noninteracting system