CH 2

2 Schrödinger equation
Our objective for the next few chapters is to learn how the Hamiltonian matrix [H] for a
given device structure (see Fig. 1.6.5) is written down. We start in this chapter with (1)
the hydrogen atom (Section 2.1) and how it led scientists to the Schrödinger equation,
(2) a simple approach called the finite difference method (Section 2.2) that can be used
to convert this differential equation into a matrix equation, and (3) a few numerical
examples (Section 2.3) showing how energy levels are calculated by diagonalizing the
resulting Hamiltonian matrix.
2.1 Hydrogen atom
Early in the twentieth century scientists were trying to build a model for atoms which
were known to consist of negative particles called electrons surrounding a positive
nucleus. A simple model pictures the electron (of mass m and charge −q) as orbiting
the nucleus (with charge Zq) at a radius r (Fig. 2.1.1) kept in place by electrostatic
attraction, in much the same way that gravitational attraction keeps the planets in orbit
around the Sun.

Zq2 mv2 Zq 2
= ⇒ v = (2.1.1)
4π ε0r 2 r 4π ε0 mr
Electrostatic force = Centripetal force
A faster electron describes an orbit with a smaller radius. The total energy of the electron
is related to the radius of its orbit by the relation
Zq 2 mv2 Zq 2
E =− + =− (2.1.2)
4π ε0r 2 8π ε0r
Potential energy + Kinetic energy = Total energy
However, it was soon realized that this simple viewpoint was inadequate since, accord-
ing to classical electrodynamics, an orbiting electron should radiate electromagnetic
waves like an antenna, lose energy continuously and spiral into the nucleus. Classically
it is impossible to come up with a stable structure for such a system except with the
33
−q
+ Zq
Fig. 2.1.1 Stationary orbits of an electron around a nucleus can be obtained by requiring their
circumferences to be integer multiples of the de Broglie wavelength.
electron sitting right on top of the nucleus, in contradiction with experiment. It was
apparent that a radical departure from classical physics was called for.
Bohr postulated that electrons could be described by stable orbits around the nucleus
at specific distances from the nucleus corresponding to specific values of angular
momenta. It was later realized that these distances could be determined by endow-
ing the electrons with a wavelike character having a de Broglie wavelength equal to
(h/mv), h being the Planck constant. One could then argue that the circumference of
an orbit had to be an integer multiple of wavelengths in order to be stable:
2πr = n(h/mv) (2.1.3)
Combining Eq. (2.1.3) with Eqs. (2.1.1) and (2.1.2) we obtain the radius and energy of
stable orbits respectively:
rn = (n 2 /Z )a0(Bohr radius) (2.1.4)

where a0 = 4π ε0 h- 2 /mq 2 = 0.0529 nm (2.1.5)
E n = −(Z /n )E 0
2 2
(2.1.6a)
where E 0 = q /8π ε0 a0 = 13.6 eV
2
(1 Rydberg) (2.1.6b)
Once the electron is in its lowest energy orbit (n = 1) it cannot lose any more energy
because there are no stationary orbits having lower energies available (Fig. 2.1.2a). If
we heat up the atom, the electron is excited to higher stationary orbits (Fig. 2.1.2b).
When it subsequently jumps down to lower energy states, it emits photons whose energy
hν corresponds to the energy difference between orbits m and n:

1 1
hν = E m − E n = E 0 Z 2
− 2 (2.1.7)
n2 m
Experimentally it had been observed that the light emitted by a hydrogen atom indeed
consisted of discrete frequencies that were described by this relation with integer values
of n and m. This striking agreement with experiment suggested that there was some
truth to this simple picture, generally known as the Bohr model.
35 2.1 Hydrogen atom
(a) Ground state (b) Excited states
n=4 n=4
n=3 n=3
n=2 n=2
n=1 n=1
Fig. 2.1.2 (a) Left to itself, the electron relaxes to its lowest energy orbit (n = 1). (b) If we heat up
the atom, the electron is excited to higher stationary orbits. When it subsequently jumps down to
lower energy states, it emits photons whose energy hν corresponds to the energy difference
between the initial and final orbits.
The Schrödinger equation put this heuristic insight on a formal quantitative basis
allowing one to calculate the energy levels for any confining potential U (
r ).
-2
∂ h
ih- = − ∇ 2 + U (
r) (2.1.8)
∂t 2m
How does this equation lead to discrete energy levels? Mathematically, one can show
r ) = −Zq 2 /4π ε0 r appropriate for a nucleus of charge
that if we assume a potential U (
+Zq, then the solutions to this equation can be labeled with three indices n, l and m
(
r , t) = φnlm ( -
r ) exp (−iE n t/h) (2.1.9)
where the energy En depends only on the index n and is given by E n = −(Z 2 /n 2 )E 0 in
agreement with the heuristic result obtained earlier (see Eq. (2.1.6a)). The Schrödinger
equation provides a formal wave equation for the electron not unlike the equation that
describes, for example, an acoustic wave in a sound box . The energy E of the electron
plays a role similar to that played by the frequency of the acoustic wave. It is well-
known that a sound box resonates at specific frequencies determined by the size and
shape of the box. Similarly an electron wave in an atomic box “resonates” at specific
energies determined by the size and shape of the box as defined by the potential energy
U (r). Let us elaborate on this point a little further.
Waves in a box: To keep things simple let us consider the vibrations u(x, t) of a
one-dimensional (1D) string described by the 1D wave equation:
∂ 2u 2∂ u
2
= v (2.1.10)
∂t 2 ∂x2
(a) (b)
x=0 x=L x=0 x=L
Fig. 2.1.3 Standing waves. (a) Acoustic waves in a “guitar” string with the displacement clamped
to zero at either end. (b) Electron waves in a one-dimensional box with the wavefunction clamped
to zero at both ends by an infinite potential.
The solutions to this equation can be written in the form of plane waves with a linear
dispersion ω = ±vk :
u = u 0 exp(ikx) exp(−iωt) ⇒ ω2 = v2 k 2 (2.1.11)
What happens if we clamp the two ends so that the displacement there is forced to be
zero (Fig. 2.1.3)? We have to superpose solutions with +k and −k to obtain standing
wave solutions. The allowed values of k are quantized leading to discrete resonant
frequencies:
u = u 0 sin(kx) exp(−iωt) ⇒ k = nπ/L ⇒ ω = nπ v/L (2.1.12)
Well, it’s the same way with the Schrödinger equation. If there is no confining potential
(U = 0), we can write the solutions to the 1D Schrödinger equation:
∂ h- 2 ∂ 2
ih- =− (2.1.13)
∂t 2m ∂ x 2
in the form of plane waves with a parabolic dispersion law E = h- 2 k2 /2m:
- ⇒ E = h- 2 k 2 /2m
= 0 exp(i kx) exp(−iEt/h) (2.1.14)
If we fix the two ends we get standing waves with quantized k and resonant frequency:
- ⇒ k = nπ/L
= 0 sin(kx) exp(−iEt/h)
⇒ E = h- 2 π 2 n 2 /2m L 2 (2.1.15)
Atomic “boxes” are of course defined by potentials U ( r ) that are more complicated
than the simple rectangular 1D potential shown in Fig. 2.1.2b, but the essential point is
the same: anytime we confine a wave to a box, the frequency or energy is discretized
because of the need for the wave to “fit” inside the box.
“Periodic” box: Another kind of box that we will often use is a ring (Fig. 2.1.4) where
the end point at x = L is connected back to the first point at x = 0 and there are no ends.
Real boxes are seldom in this form but this idealization is often used since it simplifies
the mathematics. The justification for this assumption is that if we are interested in the
properties in the interior of the box, then what we assume at the ends (or surfaces) should
37 2.1 Hydrogen atom
x = L ---
x = 0 --- x
Fig. 2.1.4 Standing waves in a ring.
make no real difference and we could assume anything that makes our calculations
simpler. However, this may not be a valid argument for “nanostructures” where the
actual surface conditions can and do affect what an experimentalist measures.
Anyway, for a periodic box the eigenfunctions are given by (cf. Eq. (2.1.15))
-
= 0 sin(kx) exp(−iEt/h)
-
and = 0 cos(kx) exp(−iEt/h)
with k = 2nπ/L ⇒ E = 2h- 2 π 2 n 2 /m L 2 (2.1.16)
The values of k are spaced by 2π/L instead of π/L, so that there are half as many
allowed values. But for each value of k there is a sine and a cosine function which have
the same eigenvalue, so that the eigenvalues now come in pairs.
An important point to note is that whenever we have degenerate eigenstates, that is,
two or more eigenfunctions with the same eigenvalue, any linear combination of these
eigenfunctions is also an eigenfunction with the same eigenvalue. So we could just as
well write the eigenstates as
-
= 0 exp(+ikx) exp(−iEt/h)
and = 0 exp(+ikx) exp(−iEt/h) -
with k = 2nπ/L ⇒ E = 2h- 2 π 2 n 2 /m L 2 (2.1.17)
This is done quite commonly in analytical calculations and the first of these is viewed
as the +k state traveling in the positive x-direction while the second is viewed as
the −k state traveling in the negative x-direction.
Electron density and probability current density: An electron with a wavefunction

(x, t) has a probability of ∗ dV of being found in a volume dV. When a number of
electrons are present we could add up ∗ for all the electrons to obtain the average
electron density n(x, t). What is the corresponding quantity we should sum to obtain
the probability current density J(x, t)?
The appropriate expression for the probability current density

ih- ∂ ∗ ∗ ∂
J= − (2.1.18)
2m ∂x ∂x
is motivated by the observation that as long as the wavefunction (x, t) obeys the
Schrödinger equation, it can be shown that
∂J ∂n
+ =0 (2.1.19)
∂x ∂t
if J is given by Eq. (2.1.18) and n = ∗ . This ensures that the continuity equation is
satisfied regardless of the detailed dynamics of the wavefunction. The electrical current
density is obtained by multiplying J by the charge (−q) of an electron.
It is straightforward to check that the “+k” and “−k” states in Eq. (2.1.17) carry
equal and opposite non-zero currents proportional to the electron density
-
J = (hk/m) ∗ (2.1.20)
-
suggesting that we associate (hk/m) with the velocity v of the electron (since we expect
J to equal nv). However, this is true only for the plane wave functions in Eq. (2.1.17).
The cosine and sine states in Eq. (2.1.16), for example, carry zero current. Indeed
Eq. (2.1.18) will predict zero current for any real wavefunction.
2.2 Method of finite differences
The Schrödinger equation for a hydrogen atom can be solved analytically, but most
other practical problems require a numerical solution. In this section I will describe one
way of obtaining a numerical solution to the Schrödinger equation. Most numerical
methods have one thing in common – they use some trick to convert the
wavefunction ( r , t) into a column vector {ψ(t)}
and the differential operator Hop into a matrix [H ]
so that the Schrödinger equation is converted from a
partial differential equation into a matrix equation
∂ d
ih- ( r , t) = Hop (
r , t) ih- {ψ(t)} = [H ] {ψ(t)}
∂t dt
This conversion can be done in many ways, but the simplest one is to choose a discrete
lattice. To see how this is done let us for simplicity consider just one dimension and
discretize the position variable x into a lattice as shown in Fig. 2.2.1: xn = na.
We can represent the wavefunction (x, t) by a column vector {ψ1 (t) ψ2 (t) · · · · · · }T
(“T” denotes transpose) containing its values around each of the lattice points at time
t. Suppressing the time variable t for clarity, we can write
{ψ1 ψ2 ··· · · ·} = {(x1 ) (x2 ) ··· · · ·}
This representation becomes exact only in the limit a → 0, but as long as a is smaller
than the spatial scale on which varies, we can expect it to be reasonably accurate.
39 2.2 Method of finite differences
x x x
x
x
x
a
Fig. 2.2.1 A continuous function can be represented by its values at a set of points on a discrete
lattice.
The next step is to obtain the matrix representing the Hamiltonian operator
h- 2 d2
Hop ≡ − + U (x)
2m dx 2
Basically what we are doing is to turn a differential equation into a difference equation.
There is a standard procedure for doing this – the finite difference technique:

∂ 2 1
→ [(xn+1 ) − 2(xn ) + (xn−1 )]
∂x2 x=xn a2
and
U (x) (x) → U (xn )(xn )
This allows us to write (note: t0 ≡ h- 2 /2ma 2 and Un ≡ U (xn ))
dψn
ih- = Hop ψ x=xn = (Un + 2t0 ) ψn − t0 ψn−1 − t0 ψn+1
dt
= [(Un + 2t0 ) δn,m − t0 δn,m+1 − t0 δn,m−1 ]ψm (2.2.1)

m
where δn,m is the Kronecker delta, which is one if n = m and zero if n = m. We can
write Eq. (2.2.1) as a matrix equation:
d
ih- {ψ(t)} = [H ] {ψ(t)} (2.2.2)
dt
The elements of the Hamiltonian matrix are given by
Hn,m = [Un + 2t0 ] δn,m − t0 δn,m+1 − t0 δn,m−1 (2.2.3)

where t0 ≡ h- 2 /2ma 2 and Un ≡ U (x n ). This means that the matrix representing H looks
like this
H= 1 2 ... N −1 N
1 2t0 + U1 −t0 0 0
2 −t0 2t0 + U2 0 0 (2.2.4)
... ... ...
N −1 0 0 2t0 + U N −1 −t0
N 0 0 −t0 2t0 + U N
For a given potential function U (x) it is straightforward to set up this matrix, once we
have chosen an appropriate lattice spacing a.
Eigenvalues and eigenvectors: Now that we have converted the Schrödinger equation
into a matrix equation (Eq. (2.2.2))
d
ih- {ψ(t)} = [H ] {ψ(t)}
dt
how do we calculate {ψ(t)} given some initial state {ψ(0)}? The standard procedure is
to find the eigenvalues E α and eigenvectors {α} of the matrix [H]:
[H ] {α} = E α {α} (2.2.5)
-
Making use of Eq. (2.2.5) it is easy to show that the wavefunction {ψ(t)} = e−iEα t/h {α}
satisfies Eq. (2.2.2). Since Eq. (2.2.2) is linear, any superposition of such solutions

-
{ψ(t)} = Cα e−iEα t/h {α} (2.2.6)
α
is also a solution. It can be shown that this form, Eq. (2.2.6), is “complete,” that is,
any solution to Eq. (2.2.2) can be written in this form. Given an initial state we can
figure out the coefficients Cα . The wavefunction at subsequent times t is then given
by Eq. (2.2.6). Later we will discuss how we can figure out the coefficients. For the
moment we are just trying to make the point that the dynamics of the system are easy
to visualize or describe in terms of the eigenvalues (which are the energy levels that we
talked about earlier) and the corresponding eigenvectors (which are the wavefunctions
associated with those levels) of [H]. That is why the first step in discussing any system
is to write down the matrix [H] and to find its eigenvalues and eigenvectors. This is
easily done using any standard mathematical package like M at l a b as we will discuss
in the next section.
41 2.3 Examples
2.3 Examples
Let us now look at a few examples to make sure we understand how to find the eigen-
energies and eigenvectors numerically using the method of finite differences described
in the last section. These examples are all simple enough to permit analytical solutions
that we can use to compare and evaluate our numerical solutions. The advantage of
the numerical procedure is that it can handle more complicated problems just as easily,
even when no analytical solutions are available.
2.3.1 Particle in a box

Consider, first the “particle in a box” problem that we mentioned in Section 2.1. The
potential is constant inside the box which is bounded by infinitely high walls at x = 0
and at x = L (Fig. 2.3.1). The eigenstates φα (x) are given by
φα (x) ∼ sin(kα x) where kα = απ/L , α = 1, 2, . . .
and their energies are given by E α = h- 2 kα2 /2m.

We could solve this problem numerically by selecting a discrete lattice with 100
points and writing down a 100 × 100 matrix [H] using Eq. (2.2.4) with all Un = 0:
H= 1 2 ... 99 100
1 2t0 −t0 0 0
2 −t0 2t0 0 0 (2.3.1)
... ... ...
99 0 0 2t0 −t0
100 0 0 −t0 2t0
U(x)
°
a = 1A
x
U=0 x x x x
1 2 99 100
Fig. 2.3.1 Energy levels for a “particle in a box” are calculated using a discrete lattice of
100 points spaced by a = 1 Å.
(b) a=1 a = 25
(a)
40 0.02
35 0.018
0.016
30
0.014
25
0.012
E (eV)
Probability
20 Analytical 0.01
15 0.008
0.006
10
0.004
5
Numerical 0.002
0 0
0 20 40 60 80 100 0 20 40 60 80 100
Eigenvalue number, a Lattice site number
Fig. 2.3.2 (a) Numerical evaluation (see Fig. 2.3.1) yields 100 eigenvalues that follow the
analytical result well for low energies but deviate at higher energies because the wavefunctions
oscillate too rapidly. (b) Probability distribution (squared eigenfunction) for eigenvalues α = 1 and
α = 25.
It is straightforward to set up this matrix and use any standard mathematical package
like M at l a b to find the eigenvalues and the corresponding eigenvectors. We obtained
100 eigenvalues, which are plotted in Fig. 2.3.2a. They follow the analytical result
E α = h- 2 π 2 α 2 /2m L 2 , with L = 101a, fairly well at low energy, but deviate at higher
energies because of the rapid oscillations in the wavefunction. Our finite difference
approximation to the second derivative operator (note that t0 ≡ h- 2 /2ma 2 )

h- 2 ∂ 2
− → t0 [(xn+1 ) − 2(xn ) + (xn−1 )]
2m ∂ x 2 x=xn
is accurate only if varies slowly enough on a length scale of a. Indeed if we put
∼ sin (kα x) it is straightforward to show that

h- 2 ∂ 2
− = t0 (kα a) 2 (xn )
2m ∂ x 2 x=xn
while
t0 [(xn+1 ) − 2(xn ) + (xn−1 )] = 2t0 (1 − cos kα a)(xn )
Since kα = απ/L, the analytical eigenvalues follow a parabolic function while the
numerical eigenvalues follow a cosine function:
E α = t0 (π αa/L)2 E α = 2t0 [1 − cos (απa/L)]
Analytical result Numerical result
The two are equivalent only if kα a = απa/L 1 so that cos (kα a) ≈ 1 − (kα2 a 2 /2).
43 2.3 Examples
Normalization: Figure 2.3.2b shows the eigenfunction squared corresponding to the

eigenvalues α = 1 and α = 25. A word about the normalization of the wavefunctions:
In analytical treatments, it is common to normalize wavefunctions such that
+∞
dx |φα (x)|2 = 1
−∞
Numerically, a normalized eigenvector satisfies the condition

N
|φα (xn )|2 = 1
n=1
So when we compare numerical results with analytical results we should expect
|φα (xn )|2 = |φα (x)|2 a (2.3.2)

Numerical Analytical
where a is the lattice constant (see Fig. 2.3.1). For example, in the present case
|φα (x)|2 = (2/L) sin2 (kα x) −→ |φα (xn )|2 = (2a/L) sin2 (kα xn )
Analytical Numerical
Since we used a = 1 Å and L = 101 Å, the numerical probability distribution should
have a peak value of 2a/L ≈ 0.02 as shown in Fig. 2.3.2b.
Boundary conditions: One more point: Strictly speaking, the matrix [H] is infinitely
large, but in practice we always truncate it to a finite number, say N, of points. This
means that at the two ends we are replacing (see Eq. (2.2.1))
−t0 ψ0 + (2t0 + U1 )ψ1 − t0 ψ2 with (2t0 + U1 )ψ1 − t0 ψ2
and
−t0 ψ N −1 + (2t0 + U N )ψ N − t0 ψ N +1 with −t0 ψ N −1 + (2t0 + U N )ψ N
In effect we are setting ψ0 and ψ N +1 equal to zero. This boundary condition is appro-
priate if the potential is infinitely large at point 0 and at point N + 1 as shown in
Fig. 2.3.3. The actual value of the potential at the end points will not affect the results
as long as the wavefunctions are essentially zero at these points anyway.
Another boundary condition that is often used is the periodic boundary condition
where we assume that the last point is connected back to the first point so that there
are no ends (Fig. 2.3.4). As we mentioned earlier (Fig. 2.1.4), the justification for this
assumption is that if we are interested in the properties in the interior of a structure,
then what we assume at the boundaries should make no real difference and we could
assume anything to make our calculations simpler.
x
0 1 2 ... N N+1
Fig. 2.3.3 The boundary condition ψ 0 = 0 and ψ N +1 = 0 can be used if we assume an infinitely
large potential at points 0 and N + 1.
x
0 1 2 ... N N+1
Fig. 2.3.4 Periodic boundary conditions assume that there are no “ends.” Point N is connected
back to point 1 as if the structure were in the form of a ring making (N + 1) equivalent to 1.
Mathematically, periodic boundary conditions are implemented by modifying the

Hamiltonian to
H= 1 2 ... 99 100
1 2t0 −t0 0 −t0
2 −t0 2t0 0 0 (2.3.3)
... ... ...
99 0 0 2t0 −t0
100 −t0 0 −t0 2t0
Note that compared to the infinite wall boundary conditions (cf. Eq. (2.3.1)) the only
change is in the elements H(1, 100) and H(100, 1). This does change the resulting
eigenvalues and eigenvectors, but the change is imperceptible if the number of points
is large. The eigenfunctions are now given by
φα (x) ∼ sin(kα x) and cos(kα x)
where kα = α 2π/L , α = 1, 2, . . . instead of
φα (x) ∼ sin(kα x)
where kα = α π/L , α = 1, 2, . . .
The values of kα are spaced by 2π /L instead of π /L, so that there are half as many
allowed values. But for each value of kα there is a sine and a cosine function which
have the same eigenvalue, so that the eigenvalues now come in pairs as evident from
Fig. 2.3.5.
45 2.3 Examples
16
14
12
10
8
E (eV)
−2
0 20 40 60 80 100
Eigenvalue number, a
Fig. 2.3.5 Energy eigenvalues for a box of length 101 Å (same as Fig. 2.3.1) with periodic
boundary conditions: the eigenvalues now come in pairs.
As we discussed earlier, instead of writing the eigenstates as
cos(kα x) and sin(kα x)
we could just as well write them as
eikα x = cos(kα x) + i sin(kα x) and e−ikα x = cos(kα x) − i sin(kα x)
This is done quite commonly in analytical calculations, but numerical calculations will
typically give the eigenvectors as cos(kα x) and sin(kα x). Both forms are equally correct
though one may be more convenient than the other for certain calculations.
Number of eigenvalues: Another important point to note about the numerical solution
is that it yields a finite number of eigenvalues (unlike the analytical solution for which
the number is infinite). This is expected since a finite matrix can have only a finite
number of eigenvalues, but one might wonder why we do not have an infinite number
of E α corresponding to an infinite number of kα a = α2πa/L , just as we have for the
analytical result. The reason is that for a discrete lattice, the wavefunctions
sin(kα x) and sin([kα + (2π/a)]x)
represent the same state because at any lattice point xn = na,
sin(kα xn ) = sin([kα + (2π/a)]xn )
They are NOT equal between two lattice points and thus represent distinct states in a
non-discrete representation. But once we adopt a discrete lattice, values of kα differing
by 2π /a represent identical states and only the values of kα a within a range of 2π yield
independent solutions. Since kα a = απa/L = απ/N , this means that there are only
N values of α that need to be considered. It is common to restrict the values of kα a to
the range (sometimes called the first Brillouin zone)
−π < kα a ≤ +π for periodic boundary conditions
and
0 < kα a ≤ +π for infinite wall boundary conditions
2.3.2 Particle in a 3D “box”

For simplicity we have limited our discussion of the method of finite differences to one
dimension, but the basic idea carries over in principle to two or three dimensions. The
diagonal elements of [H] are equal to t0 times the number of nearest neighbors (two in
one dimension, four in two dimensions and six in three dimensions) plus the potential
r ) evaluated at the lattice site, while the off-diagonal elements are equal to −t0 for
U (
neighboring sites on the lattice. That is, (ν is the number of nearest neighbors)
Hnm = νt0 n=m
= −t0 n, m are nearest neighbors (2.3.4)
=0 otherwise
However, we run into a practical difficulty in two or three dimensions. If we have lattice
points spaced by 1 Å, then a one-dimensional problem with L = 101 Å requires a matrix
[H] 100 × 100 in size. But in three dimensions this would require a matrix 106 × 106
in size. This means that in practice we are limited to very small problems. However,
if the coordinates are separable then we can deal with three separate one-dimensional
problems as opposed to one giant three-dimensional problem. This is possible if the
potential can be separated into an x-, a y-, and a z-dependent part:
r ) = Ux (x) + U y (y) + Uz (z)
U ( (2.3.5)
The wavefunction can then be written in product form:
(
r ) = X (x)Y (y)Z (z)
where each of the functions X (x), Y (y), and Z (z) is obtained by solving a separate
one-dimensional Schrödinger equation:
-2 2
h d
E x X (x) = − + Ux (x) X (x)
2m dx 2
-2 2
h d
E y Y (y) = − + U y (y) Y (y) (2.3.6)
2m dy 2
-2 2
h d
E z Z (z) = − + Uz (z) Z (z)
2m dz 2
47 2.3 Examples
The total energy E is equal to the sum of the energies associated with each of the three
dimensions: E = E x + E y + E z .
Spherically symmetric potential: Some problems may not be separable in Cartesian

coordinates but could be separable in cylindrical or spherical coordinates. For example,
the potential in a hydrogen atom U ( r ) = −q 2 /4π ε0r cannot be separated in (x, y, z).
But it is separable in (r, θ, φ) and the wavefunction may be written in the form
(r, θ, φ) = [ f (r )/r ] Ylm (θ,φ) (2.3.7)
where the radial wavefunction f (r ) is obtained by solving the radial Schrödinger

equation:
-2 2
h d l(l + 1)h- 2
E f (r ) = − + U (r ) + f (r ) (2.3.8)
2m dr 2 2mr 2
Here l = 0 for s levels, l = 1 for p levels and so on. Ylm (θ, φ) are the spherical harmonics
given by

Y00 (θ, φ) = 1/4π

Y10 (θ, φ) = 3/4π cos θ

Y1±1 (θ, φ) = ± 3/8π sin θe±iφ
etc. Equation (2.3.8) can be solved numerically using the method of finite differences
that we have described.
Normalization: Note that the overall wavefunctions are normalized such that
∞ π 2π
dr r 2
dθ sin θ dφ ||2 = 1
0 0 0
Since, from Eq. (2.3.7)
(r, θ, φ) = [f (r )/r ]Ylm (θ,φ)
and the spherical harmonics are normalized such that

π 2π
2
dθ sin θ dφ Ylm = 1
0 0
it is easy to see that the radial function f (r ) obeys the normalization condition
∞
dr |f (r )|2 = 1 (2.3.9)
0
0.06
0.025 Numerical
Numerical
0.05
0.02
Analytical
Probability
0.04
Probability
0.015
Analytical
0.03
0.01
0.02
0.005
0.01
0 0
0 1 2x (m ) --->3 4 5
0 1 2 3 4 5 x 10-10
x 10-10 r (m)
r (m)
Fig. 2.3.6 Radial probability distribution | f (r )| 2 corresponding to the two lowest eigenvalues
(−13.56 eV and −2.96 eV) for l = 0 (which correspond to the 1s and 2s levels). The dots show the
analytical result (Eqs. (2.1.10a, b)) while the solid curve denotes the numerical result obtained
using a lattice with 100 points spaced by a = 0.05 Å.
suggesting that we view |f (r )|2 as a radial probability distribution function such that
|f (r )|2 r tells us the probability of finding the electron in the volume between r and
(r +r). Numerical results with a lattice spacing of a should be compared with the
analytical values of |f (r )|2 a. For example, for the 1s and 2s levels,

|f 1s |2 a = 4ar 2 /a03 e−2r /a0 (2.3.10a)
2
r
|f 2s |2 a = ar 2 /8a03 2 − e−2r /2a0 (2.3.10b)
a0
Numerical results: If we use a lattice with 100 points spaced by a = 0.05 Å then the
two lowest eigenvalues with l = 0 (which correspond to the 1s and 2s levels) are
E 1s = −13.56 eV and E 2s = −2.96 eV
as compared with the analytical values (see Eq. (2.2.6)) E 1s = −13.59 eV and
E 2s = −3.4 eV. The 1s level agrees well, but the 2s level is considerably off. The
reason is easy to see if we plot the corresponding |f (r )|2 and compare with the analyt-
ical results. It is evident from Fig. 2.3.6 that the 1s wavefunction matches well, but it
is apparent that we do not have enough range for the 2s function. This can be fixed by
choosing a larger lattice spacing, namely a = 0.1 Å. Figure 2.3.7 shows that the wave-
function now matches the analytical result quite well and the 2s eigenvalue is −3.39 eV,
in good agreement with the analytical result. However, the 1s eigenvalue degrades
slightly to −13.47 eV, because the wavefunction is not sampled frequently enough. We
could improve the agreement for both 1s and 2s levels by using 200 points spaced by
a = 0.05 Å, so that we would have both fine sampling and large range. But the calculation
would then take longer since we would have to calculate the eigenvalues of a (200 ×
200) matrix instead of a (100 × 100) matrix.
49 Exercises
0.04
0.12
0.035
0.1
0.03
Probability
Probability
Probability --->
0.08 0.025
Probability --->
0.06 0.02
0.015
0.04
0.01
0.02
0.005
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x 10-9
x 10-9
r (m) r (m)
Fig. 2.3.7 Radial probability distribution | f (r )| 2 corresponding to the two lowest eigenvalues
(–13.47 eV and −3.39 eV) for l = 0 (which correspond to the 1s and 2s levels). Solid line shows
the analytical result (Eqs. (2.3.10a, b)) while the ×’s denote the numerical result obtained using a
lattice with 100 points spaced by a = 0.1 Å.
This simple example illustrates the essential issues one has to consider in setting up
the lattice for a numerical calculation. The lattice constant a has to be small enough to
provide adequate sampling of the wavefunction while the size of the lattice has to be
big enough to cover the entire range of the wavefunction. If it were essential to describe
all the eigenstates accurately, our problem would be a hopeless one. Luckily, however,
we usually need an accurate description of the eigenstates that lie within a certain range
of energies and it is possible to optimize our matrix [H] so as to provide an accurate
description over a desired range.
EXERCISES
E.2.1. (a) Use a discrete lattice with 100 points spaced by 1 Å to calculate the eigenen-
ergies for a particle in a box with infinite walls and compare with E α = h- 2 π 2 α 2 /2m L 2
(cf. Fig. 2.3.2a). Plot the probability distribution (eigenfunction squared) for the eigen-
values α = 1 and α = 50 (cf. Fig. 2.3.2b). (b) Find the eigenvalues using periodic
boundary conditions and compare with Fig. 2.3.5.
E.2.2. (a) Obtain the radial equation given in Eq. (2.3.8) by (1) writing the operator ∇ 2
in the Schrödinger equation in spherical coordinates:
2
∂ 2 ∂ 1 1 ∂ ∂ 1 ∂2
∇ ≡
2
+ + 2 sin θ +
∂r 2 r ∂r r sin θ ∂θ ∂θ sin2 θ ∂φ 2
(2) noting that the spherical harmonics Ylm (θ, φ) are eigenfunctions of the angular
part:

1 ∂ ∂ 1 ∂2
sin θ + Ylm = −l(l + 1)Ylm
sin θ ∂θ ∂θ sin2 θ ∂φ 2
(3) writing the wavefunction (r ) = (r ) Ylm (θ, φ) and noting that
2
∂ 2 ∂ l(l + 1)
∇ 2 = + −
∂r 2 r ∂r r2
(4) simplifying the Schrödinger equation to write for the radial part
-2 2
h ∂ 2 ∂ h- 2l(l + 1)
Eψ = − + + U (r ) + ψ
2m ∂r 2 r ∂r 2mr 2
and finally (5) writing ψ(r ) = f (r )/r , to obtain Eq. (2.3.8) for f (r ).
(b) Use a discrete lattice with 100 points spaced by a to solve Eq. (2.3.8)
-2 2
h d q2 l(l + 1)h- 2
E f (r ) = − − + f (r )
2m dr 2 4π ε0r 2mr 2
for the 1s and 2s energy levels of a hydrogen atom. Plot the corresponding radial
probability distributions |f (r )| 2 and compare with the analytical results for (a) a =
0.05 Å (cf. Fig. 2.3.6) and (b) a = 0.1 Å (cf. Fig. 2.3.7).
Strictly speaking one should replace the electron mass with the reduced mass to
account for nuclear motion, but this is a small correction compared to our level of
accuracy.
E.2.3. Use Eq. (2.1.18) to evaluate the current density associated with an electron having
the wavefunction
-
(x, t) = (e+γ x + ae−γ x )e−iEt/h
assuming γ is (a) purely imaginary (= iβ) and (b) purely real.

3 Self-consistent field
As we move from the hydrogen atom (one electron only) to multi-electron atoms, we are
immediately faced with the issue of electron–electron interactions, which is at the heart
of almost all the unsolved problems in our field. In this chapter I will explain (1) the
self-consistent field (SCF) procedure (Section 3.1), which provides an approximate way
to include electron–electron interactions into the Schrödinger equation, (2) the inter-
pretation of the energy levels obtained from this so-called “one-electron” Schrödinger
equation (Section 3.2), and (3) the energetic considerations underlying the process by
which atoms “bond” to form molecules (Section 3.3). Finally, a supplementary section
elaborates on the concepts of Section 3.2 for interested readers (Section 3.4).
3.1 The self-consistent field (SCF) procedure
One of the first successes of quantum theory after the interpretation of the hydrogen atom
was to explain the periodic table of atoms by combining the energy levels obtained from
the Schrödinger equation with the Pauli exclusion principle requiring that each level
be occupied by no more than one electron. The energy eigenvalues of the Schrödinger
equation for each value of l starting from l = 0 (see Eq. (2.3.8)) are numbered with
integer values of n starting from n = l + 1. For any (n, l) there are (2l + 1) levels with
distinct angular wavefunctions (labeled with another index m), all of which have the
same energy. For each (n, l, m) there a is an up-spin and a down-spin level making the
number of degenerate levels equal to 2(2l + 1) for a given (n, l). The energy levels look
something like Fig. 3.1.1.
The elements of the periodic table are arranged in order as the number of electrons
increases by one from one atom to the next. Their electronic structure can be written as:
hydrogen, 1s1 ; helium, 1s2 ; lithium, 1s2 2s1 ; beryllium, 1s2 2s2 ; boron, 1s2 2s2 2p1 , etc.,
where the superscript indicates the number of electrons occupying a particular orbital.
How do we calculate the energy levels for a multi-electron atom? The time-
independent Schrödinger equation
h- 2 2
E α α (
r ) = Hop α (
r ) where Hop ≡ − ∇ + U (
r)
2m
51
2 (2l + 1)
degenerate
levels
n=3
n=2
n=1
l=0 l=1 l=2
s levels p levels d levels
Fig. 3.1.1
provides a fairly accurate description of the observed spectra of all atoms, not just the
hydrogen atom. However, multi-electron atoms involve electron–electron interactions
that are included by adding a “self-consistent field (SCF),” USCF ( r ), to the nuclear
potential Unuc ( r ) = Unuc (
r ): U ( r ) + USCF (
r ), just as in Section 1.4 we added an extra
potential to the Laplace potential UL (see Eq. (1.4.1b)). The nuclear potential Unuc ,
like UL , is fixed, while USCF depends on the electronic wavefunctions and has to be
calculated from a self-consistent iterative procedure. In this chapter we will describe
this procedure and the associated conceptual issues.
Consider a helium atom consisting of two electrons bound to a nucleus with two
positive charges +2q. What will the energy levels looks like? Our first guess would be
simply to treat it just like a hydrogen atom except that the potential is
r ) = −2q 2 /4π ε0r

U (
instead of
r ) = −q 2 /4π ε0r
U (
r ) = −Zq 2 /4π ε0r we will obtain energy

If we solve the Schrödinger equation with U (
levels given by
E n = −(Z 2 /n 2 ) E 0 = −54.4 eV/n 2 (Z = 2)
just as predicted by the simple Bohr model (see Eqs. (2.1.6a, b)). However, this does not
compare well with experiment at all. For example, the ionization potential of helium
is ∼25 eV, which means that it takes a photon with an energy of at least 25 eV to ionize
a helium atom:
He + hν → He+ + e− (3.1.1a)
53 3.1 The self-consistent field (SCF) procedure
Free electron
n=3
n=2
- +2q -
n=1
Fig. 3.1.2 Ionization of a neutral helium atom takes approximately 25 eV of energy, suggesting
that the n = 1 level has an energy of −25 eV.
This suggests that the 1s level of a helium atom has an energy of −25 eV and not
−54.4 eV as the simple argument would suggest. How could we be off by over 30 eV?
It is because we did not account for the other electron in helium. If we were to measure
the energy that it takes to remove the second electron from He+
He+ + hν → He++ + e− (3.1.1b)
the result (known as the second ionization potential) is indeed close to 54.4 eV. But the
(first) ionization potential is about 30 eV less, indicating that it takes 30 eV less energy
to pull an electron out of a neutral helium atom than it takes to pull an electron out of
a helium ion (He+ ) that has already lost one electron. The reason is that an electron in
a helium atom feels a repulsive force from the other electron, which effectively raises
its energy by 30 eV and makes it easier for it to escape (Fig. 3.1.2).
In general, the ionization levels for multielectron atoms can be calculated approx-
imately from the Schrödinger equation by adding to the nuclear potential Unuc ( r ), a
self-consistent field USCF (r ) due to the other electrons (Fig. 3.1.3):
U ( r ) + USCF (
r ) = Unuc ( r) (3.1.2)
For all atoms, the nuclear potential arises from the nuclear charge of +Zq located
at the origin and is given by Unuc (r ) = −Zq 2 /4π ε0r . The self-consistent field arises
from the other (Z − 1) electrons, since an electron does not feel any potential due to
itself. In order to calculate the potential USCF (
r ) we need the electronic charge which
depends on the wavefunctions of the electron which in turn has to be calculated from
the Schrödinger equation containing USCF ( r ). This means that the calculation has to be
done self-consistently as follows.
Step 1. Guess electronic potential USCF (

r ).
Step 2. Find eigenfunctions and eigenvalues from Schrödinger equation.
Step 3. Calculate the electron density n(
r ).
Step 4. Calculate the electronic potential USCF (
r ).
(a) Nuclear charge Electronic charge
r r=0 r
r=0
Z−1
→ →
(b) Unuc(r ) USCF(r )
~ (Z−1)/r
r=0
~Z/r
r=0
Fig. 3.1.3 Sketch of (a) the nuclear charge density and the electronic charge density;
(b) potential energy felt by an additional electron due to the nucleus, Unuc (r ), and the other
electrons, USCF (r ). The latter has to be calculated self-consistently.
Step 5. If the new USCF (r ) is significantly different from last guess, update USCF (
r)
and go back to Step 2. If the new USCF ( r ) is within say 10 meV of the last
guess, the result has converged and the calculation is complete.
For Step 2 we can use essentially the same method as we used for the hydrogen atom,
although an analytical solution is usually not possible. The potential USCF ( r ) is in
general not isotropic (which means independent of θ, φ) but for atoms it can be assumed
to be isotropic without incurring any significant error. However, the dependence on r
is quite complicated so that no analytical solution is possible. Numerically, however,
it is just as easy to solve the Schrödinger equation with any U(r) as it is to solve the
hydrogen atom problem with U(r) ∼ 1/r.
For Step 3 we have to sum up the probability distributions for all the occupied
eigenstates:

f n (r ) 2
|α (
r) = r )| = r |Ylm (θ, φ)|
2 2
n( (3.1.3)
occ α occ n,l,m
55 3.1 The self-consistent field (SCF) procedure
r + r ∞
First term in Second term in

Eq. (3.1.5) Eq. (3.1.5)
Fig. 3.1.4
If we assume the charge distribution to be isotropic, we can write

σ (r ) ≡ r 2 sin θ dθ dφ n(
r) = | f n (r )|2 (3.1.4)
occ n,l,m
For Step 4 we can use straightforward electrostatics to show that

 
r 2
∞

Z −1 q 2
q σ (r )
USCF (r ) = dr σ (r ) + dr  (3.1.5)
Z 4π ε0r 4π ε0 r
0 r
The two terms in Eq. (3.1.5) arise from the contributions due to the charge within a
sphere of radius r and that due to the charge outside of this sphere as shown in Fig. 3.1.4.
The first term is the potential at r outside a sphere of charge that can be shown to be
the same as if the entire charge were concentrated at the center of the sphere:
r
q2
dr σ (r )
4π ε0r
0
The second term is the potential at r inside a sphere of charge and can be shown to
be the same as the potential at the center of the sphere (the potential is the same at all
points inside the sphere since the electric field is zero)
∞
q2 σ (r )
dr
4π ε0 r
r
We obtain the total potential by adding the two components.

To understand the reason for the factor (Z − 1)/Z in Eq. (3.1.5), we note that the
appropriate charge density for each eigenstate should exclude the eigenstate under
consideration, since no electron feels any repulsion due to itself. For example, silicon
has 14 electrons 1s2 2s2 2p6 3s2 3p2 and the self-consistent field includes all but one of
(a)
20 U SCF( r )
0
−20
U ( eV )
−40
nuc (r)
Unuc
−60
−80
−100
0 0.2 0.4 0.6 0.8 1
−9
R (m) x 10
(b)
0.1
0.09
Helium
0.08
0.07
0.06
Probability
0.05
0.04
0.03
0.02 Hydrogen
0.01
0
0 0.2 0.4 0.6 0.8 1
−9
R (m) x 10
Fig. 3.1.5 Self-consistent field method applied to the helium atom. (a) Nuclear potential Unuc (r)
and the self-consistent electronic potential USCF (r). (b) Radial probability distribution for the 1s
state in helium and hydrogen.
these electrons – for the 3p level we exclude the 3p electron, for the 3s level we exclude
the 3s electron etc. However, it is more convenient to simply take the total charge
density and scale it by the factor (Z − 1)/Z. This preserves the spherical symmetry of
the charge distribution and the difference is usually not significant. Note that the total
electronic charge is equal to Z:
∞
dr σ (r ) = 1=Z (3.1.6)
occ n,l,m
0
∞
since the radial eigenfunctions are normalized: dr |f n (r )|2 = 1.
0
Helium atom: Figure 3.1.5 shows the potential profile and the probability distribution
for the 1s state of helium obtained using the SCF method we have just described.
57 3.2 Relation to the multi-electron picture
0.08 Silicon 1s
0.07
0.06 Hydrogen 1s
0.05
Probability
0.04
0.03
0.02
Silicon 3p
0.01
0
0 1 2 3 4 5
−10
R (m) x 10
Fig. 3.1.6 Self-consistent field method applied to the silicon atom. The radial probability
distributions for hydrogen 1s level and silicon 1s level and 3p level are shown.
Also shown for comparison is the 1s level of the hydrogen atom, discussed in the
last chapter.
Silicon atom: Figure 3.1.6 shows the probability distribution for the 1s and 3p states of
silicon obtained using the SCF method. Also shown for comparison is the 1s level of the
hydrogen atom. Note that the silicon 1s state is very tightly confined relative to the 3p
state or the hydrogen 1s state. This is typical of core states and explains why such states
remain well-localized in solids, while the outer electrons (like 3p) are delocalized.
3.2 Relation to the multi-electron picture
Multi-electron Schrödinger equation: It is important to recognize that the SCF

method is really an approximation that is widely used only because the correct method
is virtually impossible to implement. For example, if we wish to calculate the eigen-
states of a helium atom with two electrons we need to solve a two-electron Schrödinger
equation of the form
-2
h
E( r1 , r2 ) = − ∇ 2 + U (
r1 ) + U (
r2 ) + Uee (
r1 , r2 ) (
r1 , r2 ) (3.2.1)
2m
where r1 and r2 are the coordinates of the two electrons and Uee is the potential energy
due to their mutual repulsion: Uee ( r1 , r2 ) = e2 /4π ε0 |
r1 − r2 |. This is more difficult to
solve than the “one-electron” Schrödinger equation that we have been talking about,
but it is not impossible. However, this approach quickly gets out of hand as we go to
bigger atoms with many electrons and so is seldom implemented directly. But suppose
we could actually calculate the energy levels of multi-electron atoms. How would we
use our results (in principle, if not in practice) to construct a one-electron energy level
Energy levels of Energy levels of

(N − 1)-electron atom Lowest energy (N + 1)-electron atom
level of N- electron atom
E n(N − 1) E n(N + 1)
Ionization E G (N) Affinity
Fig. 3.2.1 One-electron energy levels represent energy differences between the energy levels of
the N-electron atom and the (N − 1)- or the (N + 1)-electron atom. The former (called the
ionization levels) are the filled states from which an electron can be removed while the latter
(the affinity levels) are the empty states to which an electron can be added.
diagram like the ones we have been drawing? The answer depends on what we want
our one-electron energy levels to tell us.
Ionization levels and affinity levels: Our interest is primarily in describing the flow of
current, which involves inserting an electron and then taking it out or vice versa, as we
discussed in Chapter 1. So we would want the one-electron energy levels to represent
either the energies needed to take an electron out of the atom (ionization levels) or the
energies involved in inserting an electron into the atom (affinity levels) (Fig. 3.2.1).
For the ionization levels, the one-electron energies εn represent the difference
between the ground state energy E G (N ) of the neutral N-electron atom and the
nth energy level E n (N − 1) of the positively ionized (N − 1)-electron atom:
εn = E G (N ) − E n (N − 1) (3.2.2a)
These ionization energy levels are measured by looking at the photon energy needed to
ionize an electron in a particular level. Such photoemission experiments are very useful
for probing the occupied energy levels of atoms, molecules, and solids. However, they
only provide information about the occupied levels, like the 1s level of a helium atom
or the valence band of a semiconductor. To probe unoccupied levels such as the 2s
level of a helium atom or the conduction band of a semiconductor we need an inverse
photoemission (IPE) experiment (see Fig. 3.2.2):
He + e− → He− + hν
with which to measure the affinity of the atom for acquiring additional electrons. To
calculate the affinity levels we should look at the difference between the ground state
energy E G (N ) and the nth energy level E n (N + 1) of the negatively ionized (N + 1)-
electron atom:
εn = E n (N + 1) − E G (N ) (3.2.2b)
Vacuum
level
Electron affinity, EA
Ionization potential, IP
Ionization Affinity
levels levels
Fig. 3.2.2 The ionization levels include the repulsive potential from Z − 1 electrons while the
affinity levels include that of Z electrons, so that the latter is higher in energy by the single-electron
charging energy U0 .
Note that if we want the energy levels to correspond to optical transitions then we should
look at the difference between the ground state energy E G (N ) and the nth energy level
E n (N ) of the N-electron atom, since visible light does not change the total number of
electrons in the atom, just excites them to a higher energy:
εn = E n (N ) − E G (N )
There is no a priori reason why the energy gap obtained from this calculation should
correspond to the energy gap obtained from either the ionization or the affinity levels.
In large solids (without significant excitonic effects) we are accustomed to assuming
that the optical gap is equal to the gap between the valence and conduction bands, but
this need not be true for small nanostructures.
Single-electron charging energy: As we have explained above, the straightforward

approach for calculating the energy levels would be to calculate the energies E G (N )
and E n (N ± 1) from an N-electron and an (N ± 1)-electron Schrödinger equation
(cf. Eq. (3.2.1) which is a two-electron Schrödinger equation) respectively. This, how-
ever, is usually impossible and the only practical approach for large atoms, molecules,
or solids is to include an effective potential USCF (
r ) in the Schrödinger equation as we
have been discussing.
How do we choose this effective potential? If we use Uee (N ) to denote the total
electron–electron interaction energy of an N-electron system then the appropriate USCF
for the ionization levels is equal to the change in the interaction energy as we go from
an N-electron to an (N − 1)-electron atom:
[USCF ]ionization = Uee (N ) − Uee (N − 1) (3.2.3a)
Similarly the appropriate USCF for the affinity levels is equal to the change in the
interaction energy between an N-electron and an (N + 1)-electron atom:
[USCF ]affinity = Uee (N + 1) − Uee (N ) (3.2.3b)
The electron–electron interaction energy of a collection of N electrons is proportional

to the number of distinct pairs:
Uee (N ) = U0 N (N − 1)/2 (3.2.4)
where U0 is the average interaction energy per pair, similar to the single-electron charg-
ing energy introduced in Section 1.4. From Eqs. (3.2.3a, b) and (3.2.4) it is easy to see
that
[USCF ]ionization = U0 (N − 1) while [USCF ]affinity = U0 N (3.2.5)
This means that to calculate the ionization levels of a Z-electron atom, we should use
the potential due to (Z − 1) electrons (one electron for helium) as we did in the last
section. But to calculate the affinity levels we should use the potential due to Z electrons
(two electrons for helium). The energy levels we obtain from the first calculation are
lower in energy than those obtained from the second calculation by the single-electron
charging energy U0 .
As we discussed in Section 1.5, the single-electron charging energy U0 depends on
the degree of localization of the electronic wavefunction and can be several electron-
volts in atoms. Even in nanostructures that are say 10 nm or less in dimension, it can
be quite significant (that is, comparable to kB T).
Typically one uses a single self-consistent potential
USCF = ∂Uee /∂ N = U0 N − (U0 /2) (3.2.6)
for all levels so that the ionization levels are (U0 /2) lower while the affinity levels are
(U0 /2) higher than the energy levels we calculate. One important consequence of this is
that even if an SCF calculation gives energy levels that are very closely spaced compared
to kB T (see Fig. 3.2.3a), a structure may not conduct well, because the one-electron
charging effects will create a “Coulomb gap” between the occupied and unoccupied
levels (Fig. 3.2.3b). Of course, this is a significant effect only if the single-electron
charging energy U0 is larger than kB T.
Hartree approximation: In large conductors (large R) U0 is negligible and the dis-

tinction between Z and (Z − 1) can be ignored. The self-consistent potential for both
ionization and affinity levels is essentially the same and the expression
USCF = ∂Uee /∂ N (3.2.7)

(a) (b)
Fig. 3.2.3
can be generalized to obtain the standard expression used in density functional theory
(DFT):
∂Uee
r) =
USCF ( (3.2.8)
∂[n(
r )]
which tells us that the self-consistent potential at any point r is equal to the change in
the electron–electron interaction energy due to an infinitesimal change in the number
of electrons at the same point. If we use the standard expression for Uee from classical
electrostatics

1 q 2 n( r )
r )n(
Uee = dr d r (3.2.9)
2 4π ε | r − r |
Equation (3.2.8) yields the Hartree approximation, UH (
r ) for the self-consistent
potential:

r )
q 2 n(
r ) = d
UH ( r (3.2.10)
r − r |
4π ε |
which is a solution of the Poisson equation −∇ 2 UH = −q 2 n/ε in a homogeneous
medium. Device problems often require us to incorporate complicated boundary con-
ditions including different materials with different dielectric constants. It is then more
convenient to solve a modified form of the Poisson equation that allows a spatially
varying relative permittivity:
· (εr∇UH) = q 2 n/ε0
−∇ (3.2.11)
But for atoms, there is no complicated inhomogeneity to account for and it is more
convenient to work with Eq. (3.2.10).
Correlation energy: The actual interaction energy is less than that predicted by
Eq. (3.2.9) because electrons can correlate their motion so as to avoid each other –
this correlation would be included in a many-electron picture but is missed in the
one-particle picture. One way to include it is to write

1 e2 n( r ) [1 − g(
r )n( r , r )]
Uee = d
r d r
2 4π ε | r − r |
where g is a correlation function that accounts for the fact that the probability of
finding two electrons simultaneously at r and r is not just proportional to n( r ), but
r )n(
is somewhat reduced because electrons try to avoid each other (actually this correlation
factor is spin-dependent, but we are ignoring such details). The corresponding self-
consistent potential is also reduced (cf. Eq. (3.2.10)):

r ) [1 − g(
e2 n( r , r )]
USCF = d r (3.2.12)
4π ε | r − r |
Much research has gone into estimating the function g( r , r ) (generally referred to as
the exchange-correlation “hole”).
The basic effect of the correlation energy is to add a negative term Uxc ( r ) to the
Hartree term UH (
r ) discussed above (cf. Eq. (3.2.10)):
r ) = UH (
USCF ( r ) + Uxc (
r) (3.2.13)
One simple approximation, called the local density approximation (LDA) expresses
Uxc at a point in terms of the electron density at that point:
q2
r) = −
Uxc ( r )] 1/3
C[n( (3.2.14)
4π ε0
Here, C is a constant of order one. The physical basis for this approximation is that
an individual electron introduced into a medium with a background electron density
n(r) will push other electrons in its neighborhood, creating a positive correlation “hole”
around it. If we model this hole as a positive sphere of radius r0 then we can estimate
r0 by requiring that the total charge within the sphere be equal in magnitude to that of
an electron:
1
n(r ) 4πr03 /3 = 1 → r0 = [n(r )]−1/3
C
C being a constant of order one. The potential in Eq. (3.2.14) can be viewed as the
potential at the center of this positive charge contained in a sphere of radius r0 :
q2
r) = −
Uxc (
4π ε0 r0
Much work has gone into the SCF theory and many sophisticated versions of Eq. (3.2.14)
have been developed over the years. But it is really quite surprising that the one-electron
picture with a suitable SCF often provides a reasonably accurate description of multi-
electron systems. The fact that it works so well is not something that can be proved
mathematically in any convincing way. Our confidence in the SCF method stems from
the excellent agreement that has been obtained with experiment for virtually every atom
in the periodic table (see Fig. 3.2.4). Almost all the work on the theory of electronic
structure of atoms, molecules, and solids is based on this method and that is what we
will be using.
63 3.3 Bonding
−10−1 −10−1
Experimental Experimental
Calculated Calculated
7s
−100 −100
6s
6p
5s
−101 −101 5p
E (Rydbergs)
4s E (Rydbergs) 4p
−102 −102
3s 3p
2s 2p
−103 −103
1s
−104 −104
0 20 40 60 80 100 0 20 40 60 80 100
Atomic number, Z Atomic number, Z
Fig. 3.2.4 Energy levels as a function of the atomic number calculated theoretically using a
self-consistent field method. The results are in excellent agreement with experiment (adapted from
Herman and Skillman (1963)). For a hydrogen atom, the s and p levels are degenerate (that is, they
have the same energy). This is a consequence of the ∼1 / r dependence of the nuclear potential. But
this is not true of the self-consistent potential due to the electrons and, for multi-electron atoms, the
s state has a lower energy than the p state.
3.3 Bonding
One of the first successes of quantum theory was to explain the structure of the periodic
table of atoms by combining the energy levels obtained from the Schrödinger equation
with the Pauli exclusion principle requiring that each level be occupied by no more
than one electron. In Section 3.3.1 we will discuss the general trends, especially the
periodic character of the energy levels of individual atoms. We will then discuss two
bonding mechanisms (ionic (Section 3.3.2) and covalent (Section 3.3.3)) whereby a
pair of atoms, A and B, can lower their overall energy by forming a molecule AB:
E(AB) < E(A) + E(B).
−10
Energy (eV)
−20
−30
−40
−50
0 20 40 60 80 100
Atomic number
Fig. 3.3.1 Energy of the outermost s (×) and p levels (◦) of the first 86 elements of the periodic
table excluding the d- and f-shell transition metals (Z = 21–28, 39–46, and 57–78). The numbers
are taken from Harrison (1999) and Mann (1967).
3.3.1 Valence electrons

It is important to note that only the electrons in the outermost shell, referred to as the
valence electrons, participate in the bonding process. The energies of these valence
electrons exhibit a periodic variation as shown in Fig. 3.3.1 for the first 86 atoms of
the periodic table from hydrogen (atomic number Z = 1) to radon (Z = 86), exclud-
ing the d- and f-shell transition metals (see Table 3.3.1). The main point to notice is
that the energies tend to go down as we go across a row of the periodic table from
lithium (Li) to neon (Ne), increase abruptly as we step into the next row with sodium
(Na) and then decrease as we go down the row to argon (Ar). This trend is shown by
both the s and p levels and continues onto the higher rows. Indeed this periodic variation
in the energy levels is at the heart of the periodic table of the elements.
3.3.2 Ionic bonds

Ionic bonds are typically formed between an atom to the left of the periodic table (like
sodium, Na) and one on the right of the periodic table (like chlorine, Cl). The energy
levels of Na and Cl look roughly as shown in Fig. 3.3.2. It seems natural for the 3s
electron from Na to “spill over” into the 3p levels of Cl, thereby lowering the overall
energy as shown. Indeed it seems “obvious” that the binding energy, Ebin , of NaCl
would be
E bin = E(Na) + E(Cl) − E(Na+ Cl− ) = 12.3 − 5.1 = 7.2 eV.

Table 3.3.1 First 86 atoms of the periodic table from hydrogen (atomic number Z = 1) to radon (Z = 86), excluding the d- and f-shell
transition metals (Z = 21–28, 39–46, and 57–78)
H He
(Z=1) (Z=2)
Li Be B C N O F Ne
(Z=3) (Z=4) (Z=5) (Z=6) (Z=7) (Z=8) (Z=9) (Z=10)
Na Mg Al Si P S Cl Ar
(Z=11) (Z=12) (Z=13) (Z=14) (Z=15) (Z=16) (Z=17) (Z=18)
K Ca ... Cu Zn Ga Ge As Se Br Kr
(Z=19) (Z=20) ... (Z=29) (Z=30) (Z=31) (Z=32) (Z=33) (Z=34) (Z=35) (Z=36)
Rb Sr ... Ag Zn In Sn Sb Te I Xe
(Z=37) (Z=38) ... (Z=47) (Z=48) (Z=49) (Z=50) (Z=51) (Z=52) (Z=53) (Z=54)
Cs Ba ... Au Hg Tl Pb B Po At Rn
(Z=55) (Z=56) ... (Z=79) (Z=80) (Z=81) (Z=82) (Z=83) (Z=84) (Z=85) (Z=86)
Na Cl + -
Na Cl
3s 3s
−5 eV −5 eV
3p 3p
−13.8 eV −13.8 eV
3s 3s
−29.2 eV −29.2 eV
Fig. 3.3.2 Formation of Na+ Cl− from individual Na and Cl atoms with a 3s electron from Na
“spilling over” into the 3p levels of Cl thereby lowering the overall energy. This is only part of the
story, since the overall energetics also includes the electrostatic energy stored in the microscopic
capacitor formed by the two ions as explained in the text.
But this argument is incomplete because we also need to consider the change in
the electrostatic energy due to the bonding. The correct binding energy is more
like 4 eV.
The point is that the energy levels we have drawn here are all ionization levels. The
energy needed to create a sodium ion is given by its ionization potential (IP)
E(Na+ ) − E(Na) = IP(Na) = 5 eV (3.3.1a)
But the energy needed to create a chlorine ion is given by the electron affinity (EA) of
Cl and this includes an extra charging energy U0 :
E(Cl) − E(Cl− ) = EA(Cl) = IP(Cl) − U0 = 13.8 eV − U0 (3.3.1b)
Combining Eqs. (3.3.1a) and (3.3.1b) we obtain
E(Na) + E(Cl) − E(Na+ ) − E(Cl− ) = 8.8 eV − U0 (3.3.2)
However, this is not the binding energy of NaCl. It gives us the energy gained in
converting neutral Na and neutral Cl into a Na+ and a Cl− ion completely separated
from each other. If we let a Na+ and a Cl− ion that are infinitely far apart come together
to form a sodium chloride molecule, Na+ Cl− , it will gain an energy U0 in the process.
E(Na+ ) + E(Cl− ) − E(Na+ Cl− ) = U0
so that the binding energy is given by
E bin = E(Na) + E(Cl) − E(Na+ Cl− ) = 8.8 eV − U0 + U0 (3.3.3)
U0 − U0 is approximately 5 eV, giving a binding energy of around 4 eV. The numerical
details of this specific problem are not particularly important or even accurate. The
main point I wish to make is that although the process of bonding by electron transfer
may seem like a simple one where one electron “drops” off an atom into another with
67 3.3 Bonding
H H
H H
EA
1s 1s
E0 E0
E0= −13.6 eV EB
Fig. 3.3.3 Formation of H2 from individual H atoms with a bonding level EB and an anti-bonding
level EA .
a lower energy level, the detailed energetics of the process require a more careful
discussion. In general, care is needed when using one-electron energy level diagrams
to discuss electron transfer on an atomic scale.
3.3.3 Covalent bonds

We have just seen how a lowering of energy comes about when we bring together
an atom from the left of the periodic table (like sodium) and one from the right (like
chlorine). The atoms on the right of the periodic table have lower electronic energy
levels and are said to be more electronegative than those on the left. We would expect
electrons to transfer from the higher energy levels in the former to the lower energy
levels in the latter to form an ionic bond.
However, this argument does not explain covalent bonds which involve atoms with
roughly the same electronegativity. The process is a little more subtle. For example, it
is hard to see why two identical hydrogen atoms would want to form a H2 molecule,
since no lowering of energy is achieved by transferring an electron from one atom to
the other. What happens is that when the two atoms come close together the resulting
energy levels split into a bonding level (EB ) and an anti-bonding level (EA ) as shown
in Fig. 3.3.3. Both electrons occupy the bonding level which has an energy lower than
that of an isolated hydrogen atom: EB < E0 .
How do we calculate EB ? By solving the Schrödinger equation:
-2
h
E α α (
r) = − ∇ + UN (
2
r ) + UN (
r ) + USCF (
r ) α (
r) (3.3.4)
2m
where UN (r ) and UN (r ) are the potentials due to the left and the right nuclei respectively
and USCF (r) is the potential that one electron feels due to the other. To keep things simple
let us ignore USCF (r) and calculate the electronic energy levels due to the nuclear
potentials alone:
-2
h
E α0 α0 (
r) = − ∇ + UN (
2
r ) + UN (
r ) α0 (r) (3.3.5)
2m
25
U N ,N '
20
15
10
Ue, e'
5
Energy (eV)
−5
E B0 − E 0
−10
−15
−20 R0
−25
0 1 2 3 4
R (m) −10
x 10
Fig. 3.3.4 Various energies as a function of the nuclear distance R. ×××, approximate
electron–electron repulsive energy (Ue,e ). Solid curve, nucleus–nucleus repulsive energy (UN, N ).
Dashed curve, EB0 −E0 ; energy of the bonding level in a H2 molecule relative to the 1s level in a
hydrogen atom calculated approximately from the Schrödinger equation without any
self-consistent potential. ++++, binding energy of a H2 molecule relative to two hydrogen atoms
estimated from 2(EB0 −E0 ) + UN,N + Ue,e .
The lowest energy solution to Eq. (3.3.5) can be written approximately as

a+b
E B0 = E 0 + (3.3.6)
1+s
where
1 − (1 + R)e−2R
a = −2E 0
R
b = −2E 0 (1 + R) e−R
2
s = e−R [1 + R + (R /3)]
R ≡ R/a0
R being the center-to-center distance between the hydrogen atoms.

Let us now try to understand the competing forces that lead to covalent bonding.
The dashed curve in Fig. 3.3.4 shows EB0 − E0 versus the bond length R as given
by Eq. (3.3.6). Experimentally, the bond length R for a H2 molecule is 0.074 nm,
indicating that the overall energy is a minimum for this value of R. Since the energy
keeps decreasing as R is decreased, one might wonder why the two hydrogen atoms do
not just sit on top of each other (R = 0). To answer this question we need to calculate the
overall energy which should include the electron–electron repulsion (note that USCF (r)
69 3.3 Bonding
was left out from Eq. (3.3.6)) as well as the nucleus–nucleus repulsion. To understand
the overall energetics let us consider the difference in energy between a hydrogen
molecule (H2 ) and two isolated hydrogen atoms (2H).
The energy required to assemble two separate hydrogen atoms from two protons
(N, N ) and two electrons (e, e ) can be written as
E(2H) = Ue,N + Ue ,N = 2E 0 (3.3.7a)
The energy required to assemble an H2 molecule from two protons (N, N ) and two
electrons (e, e ) can be written as
E(H2 ) = UN,N + Ue,e + Ue,N + Ue,N + Ue ,N + Ue ,N (3.3.7b)
Equation (3.3.6) gives the quantum mechanical value of (Ue,N + Ue,N ) as well as
(Ue ,N + Ue ,N ) as E B0 . Hence
E(H2 ) = UN,N + Ue,e + 2E B0 (3.3.7c)
The binding energy is the energy it takes to make the hydrogen molecule dissociate
into two hydrogen atoms and can be written as
E bin = E(H2 ) − E(2H) = 2(E B0 − E 0 ) + UN,N + Ue,e (3.3.8)
This is the quantity that ought to be a minimum at equilibrium and it consists of three
separate terms. Eq. (3.3.6) gives us only the first term. The second term is easily written
down since it is the electrostatic energy between the two nuclei, which are point charges:
UN,N = q 2 /4π ε0 R (3.3.9a)
The electrostatic interaction between the two electrons should also look like q 2 /4π ε0 R
for large R, but should saturate to ∼ q 2 /4π ε0 a0 at short distances since the electronic
charges are diffused over distances ∼a0 . Let us approximate it as

Ue,e ∼
= q 2 /4π ε0 R 2 + a02 (3.3.9b)
noting that this is just an oversimplified approximation to what is in general a very dif-
ficult quantum mechanical problem – indeed, electron–electron interactions represent
the central outstanding problem in the quantum theory of matter.
The solid curve in Fig. 3.3.4 shows UN,N (Eq. (3.3.9a)), while the ××× curve shows
Ue,e (Eq. (3.3.9b)). The +++ curve shows the total binding energy estimated from
Eq. (3.3.8). It has a minimum around 0.1 nm, which is not too far from the experimental
bond length of 0.074 nm. Also the binding energy at this minimum is ∼4.5 eV, very
close to the actual experimental value. Despite the crudeness of the approximations
used, the basic physics of bonding is illustrated fairly well by this example.
H H
Fig. 3.3.5 A hydrogen molecule can be viewed as two masses connected by a spring.
Vibrational frequency: The shape of the binding energy vs. R curve suggests that
we can visualize a hydrogen molecule as two masses connected by a spring (Fig.
3.3.5). An ideal spring with a spring constant K has a potential energy of the form
U (R) = K (R − R0 )2 /2. The binding energy of the hydrogen molecule (see Fig. 3.3.4)
can be approximated as U (R) ∼= U (R0 ) + K (R − R0 )2 /2, where the effective spring
constant K is estimated from the curvature [d2 U/dR 2 ] R=R0 . Indeed the vibrational
√
frequency of the H–H bond can be estimated well from the resonant frequency 2K /M
of the mass and spring system where M is the mass of a hydrogen atom.
Ionization levels: As we have discussed, the energy levels of a multi-electron system

usually denote the ionization levels, that is the energy it takes to strip an electron from
the system. This means that in the present context the energy level EB for a hydrogen
molecule should represent
E B = E(H2 ) − E(H+
2)
Since E(H+
2 ) = UN,N + Ue ,N + Ue ,N , we can write using Eq. (3.3.7b),

E B = Ue,e + Ue,N + Ue,N = Ue,e + E B0 (3.3.10)
It is easy to check that for our model calculation (see Fig. 3.3.4) EB0 is nearly 15 eV
below E0 , but EB lies only about 4 eV below E0 . If we were to include a self-consistent
field USCF (r) in the Schrödinger equation, we would obtain the energy EB which would
be higher (less negative) than the non-interacting value of EB0 by the electron–electron
interaction energy Ue,e .
Binding energy: It is tempting to think that the binding energy is given by
E bin = 2(E B − E 0 ) + UN,N
since EB includes the electron–electron interaction energy Ue,e . However, it is easy to

see from Eqs. (3.3.8) and (3.3.10) that the correct expression is
E bin = 2(E B − E 0 ) + (UN,N − Ue,e )
The point I am trying to make is that if we include the electron–electron interaction

in our calculation of the energy level EB then the overall energy of two electrons
is NOT 2EB , for that would double-count the interaction energy between the two
71 3.4 Supplementary notes: multi-electron picture
electrons. The correct energy is obtained by subtracting off this double-counted part:
2EB − Ue,e .
3.4 Supplementary notes: multi-electron picture
As I mentioned in Section 3.2, the SCF method is widely used because the exact method
based on a multi-electron picture is usually impossible to implement. However, it is
possible to solve the multi-electron problem exactly if we are dealing with a small
channel weakly coupled to its surroundings, like the one-level system discussed in
Section 1.4. It is instructive to recalculate this one-level problem in the multi-electron
picture and compare with the results obtained from the SCF method.
One-electron vs. multi-electron energy levels: If we have one spin-degenerate level

with energy ε, the one-electron and multi-electron energy levels would look as shown
in Fig. 3.4.1. Since each one-electron energy level can either be empty (0) or occupied
(1), multi-electron states can be labeled in the form of binary numbers with a number
of digits equal to the number of one-particle states. N one-electron states thus give
rise to 2N multi-electron states, which quickly diverges as N increases, making a direct
treatment impractical. That is why SCF methods are so widely used, even though they
are only approximate.
Consider a system with two degenerate one-electron states (up-spin and down-spin)
that can either be filled or empty. All other one-electron states are assumed not to
change their occupation: those below remain filled while those above remain empty.
Let us assume that the electron–electron interaction energy is given by
Uee (N ) = (U0 /2)N (N − 1) (same as Eq. (3.2.4))
One-electron energy levels Multi-electron energy levels
11
m1 E0 + e + (U0 /2)
e
m2
01 10
Source Drain E0
V
E0 − e + (U0/2)
I I 00
Fig. 3.4.1 One-electron vs. multi-electron energy levels in a channel with one spin-degenerate
level having energy ε.
corresponding to a self-consistent potential (see Eq. (3.2.6))
∂Uee /∂ N = U0 N − (U0 /2)
Suppose the number of electrons N0 in the neutral state corresponds to having one of
these states filled. The one-electron energy levels ε can be written as the sum of the
“bare” levels ε̃ (obtained from a Schrödinger equation with just the nuclear potential,
UN ) plus the self-consistent potential [∂Uee /∂ N ] N =N0 :
ε = ε̃ + [∂Uee /∂ N ] N =N0 = ε̃ + U0 N0 − (U0 /2)
Consider now the multi-electron picture. We have four available multi-electron states
which we can designate as 00, 01, 10, and 11. In the neutral state, the system is in either
the (10) or the (01) state whose total energy we denote as
E(10) = E(01) ≡ E 0
We can write the energies of the other multi-electron states as
E(11) = E 0 + ε̃ + Uee (N0 + 1) − Uee (N0 )

= E 0 + ε̃ + U0 N0 = E 0 + ε + (U0 /2)
and
E(00) = E 0 − ε̃ − Uee (N0 ) + Uee (N0 − 1)

= E 0 − ε̃ − U0 (N0 − 1) = E 0 − ε + (U0 /2)
Master equation: In the multi-electron picture, the overall system has different proba-
bilities Pα of being in one of the 2N possible states α and all the probabilities must add
up to one:

Pα = 1 → P00 + P01 + P10 + P11 = 1 (3.4.1)

α
We can calculate the individual probabilities by noting that the system is continually
shuffled among these states and under steady-state conditions there must be no net flow
into or out of any state:

R (α → β) Pα = R (β → α)Pβ (3.4.2)
β β
Knowing the rate constants, we can calculate the probabilities by solving Eq. (3.4.2).
Equations involving probabilities of different states are called master equations. We
could call Eq. (3.4.2) a multi-electron master equation.
The rate constants R(α → β) can be written down assuming a specific model for the
interaction with the surroundings. For example, if we assume that the interaction only
involves the entry and exit of individual electrons from the source and drain contacts
then for the 00 and 01 states the rate constants are given by
01 E0 
γ1 γ2  γ γ2
  1
f 1 + f 2   (1 − f 1 ) + (1 − f 2 )
h h  h h
E 0 − ε + (U0 /2)
00
where
f 1 ≡ f 0 (ε1 − µ1 ) and f 2 ≡ f 0 (ε1 − µ2 )
tell us the availability of electrons with energy ε1 = ε − (U0 /2) in the source and drain
contacts respectively. The entry rate is proportional to the available electrons, while
the exit rate is proportional to the available empty states. The same picture applies
to the flow between the 00 and the 10 states, assuming that up- and down-spin states
are described by the same Fermi function in the contacts, as we would expect if each
contact is locally in equilibrium.
Similarly we can write the rate constants for the flow between the 01 and the 11
states
11 E 0 − ε + (U0 /2) 
γ1 γ2  γ γ2
  1
f 1 + f 2   (1 − f 1 ) + (1 − f 2 )
h h  h h
E0
01
where
f 1 ≡ f 0 (ε2 − µ1 ) and f 2 ≡ f 0 (ε2 − µ2 )
tell us the availability of electrons with energy ε2 = ε + (U0 /2) in the source and drain
contacts corresponding to the energy difference between the 01 and 11 states. This is
larger than the energy difference ε between the 00 and 01 states because it takes more
energy to add an electron when one electron is already present due to the interaction
energy U0 .
Using these rate constants it is straightforward to show from Eq. (3.4.2) that
P10 P01 γ1 f 1 + γ2 f 2
= = (3.4.3a)
P00 P00 γ1 (1 − f 1 ) + γ2 (1 − f 2 )
and
P11 P11 γ1 f 1 + γ2 f 2
= = (3.4.3b)
P10 P01 γ1 (1 − f 1 ) + γ2 (1 − f 2 )
Together with Eq. (3.4.1), this gives us all the individual probabilities. Figure 3.4.2
shows the evolution of these probabilities as the gate voltage VG is increased
Gate, VG
m e + UL m
Source Drain
e2 + UL Evolution of
e1 + U L one-electron
energy levels
m
e1 ≡ e − (U0 /2)
e1 + U L
e 2 + UL e2 ≡ e + (U0 /2)
00 11
1
0.9
0.8
Evolution of
Probability
0.7
0.6
many-electron
10
0.5 state with gate
0.4 01 voltage
0.3
0.2
0.1
0
0 0.2 0.4 0.6 0.8 1
Gate voltage, VG (V)
Fig. 3.4.2 Evolution of the energy levels of a channel with one spin-degenerate level as the gate
voltage VG is made more positive, holding the drain voltage VD equal to zero. µ = 0, ε = 0.2 eV,
kB T = 0.025 eV, U0 = 0.25 eV, UL = −qVG . Lower plot shows the probabilities of finding the
channel in one of its four states: P00 (°), P01 = P10 (solid) and P11 (×).
holding the drain voltage VD equal to zero. The gate voltage shifts the one-electron level
ε → ε + UL (we have assumed UL = −qVG ) and the probabilities are calculated from
Eqs. (3.4.3a, b) and (3.4.1) noting that the Fermi functions are given by
f 1 = f 0 (ε1 + UL − µ1 ), f 2 = f 0 (ε1 + UL − µ2 ) (3.4.4a)

f 1 = f 0 (ε2 + UL − µ1 ), f 2 = f 0 (ε2 + UL − µ2 ) (3.4.4b)
The system starts out in the 00 state (P00 = 1), shifts to the 01 and 10 states (P01 =
P10 = 0.5) once ε1 + UL drops below µ, and finally goes into the 11 state (P11 = 1)
when ε2 + UL drops below µ.
Relation between the multi-electron picture and the one-electron levels: As I have
emphasized in Section 3.2, one-electron energy levels represent differences between
energy levels in the multi-electron picture corresponding to states that differ by one
electron. Transitions involving the addition of one electron are called affinity levels
while those corresponding to the removal of one electron are called ionization levels.
For example (see Fig. 3.4.2), if the system is in the 00 state then there are two degenerate
one-electron levels ε1 + UL representing
ε1 + U L = E(10) − E(00) = E(01) − E(00) Affinity levels
Once it is in the 10 state there are two one-electron levels
ε1 + UL = E(10) − E(00) Ionization level

and ε2 + UL = E(11) − E(10) Affinity level
In the 11 state there are two degenerate one-electron levels
ε2 + UL = E(11) − E(10) = E(11) − E(01) Ionization levels
Affinity levels lie above µ, while ionization levels lie below µ as shown in Fig. 3.4.2.
This is a very important general concept regarding the interpretation of the one-electron
energy levels when dealing with complicated interacting objects. The occupied (or
ionization) levels tell us the energy levels for removing an electron while the unoccupied
(or affinity) levels tell us the energy levels for adding an extra electron. Indeed that is
exactly how these levels are measured experimentally, the occupied levels by photo-
emission (PE) and the unoccupied levels by inverse photoemission (IPE) as mentioned
in Section 1.1.
Law of equilibrium: Figure 3.4.2 represents an equilibrium calculation with both

source and drain contacts having the same Fermi function: f1 = f2 . Equilibrium problems
do not really require the use of a master equation like Eq. (3.4.2). We can use the general
principle of equilibrium statistical mechanics which states that the probability Pα that
the system is in a multi-electron state α with energy Eα and Nα electrons is given by
1
Pα = exp[−(E α − µNα )/kB T ] (3.4.5)
Z
where the constant Z (called the partition function) is determined so as to ensure that
the probabilities given by Eq. (3.4.5) for all states α add up to one:

Z= exp[−(E α − µNα )/kB T ] (3.4.6)

α
This is the central law of equilibrium statistical mechanics that is applicable to any
system of particles (electrons, photons, atoms, etc.), interacting or otherwise (see for
example, Chapter 1 of Feynman, 1972). The Fermi function is just a special case of
this general relation that can be obtained by applying it to a system with just a single
one-electron energy level, corresponding to two multi-electron states:
α Nα Eα Pα
0 0 0 1/Z
1 1 ε (1/Z ) exp[(µ − ε)/kB T ]
so that Z = 1 + exp[(µ − ε)/kB T ] and it is straightforward to show that the average
number of electrons is equal to the Fermi function (Eq. (1.1.1)):

exp[(µ − ε)/kB T ]
N= Nα Pα = P1 = = f 0 [ε − µ]
α 1 + exp[(µ − ε)/kB T ]
For multi-electron systems, we can use the Fermi function only if the electrons are
not interacting. It is then justifiable to single out one level and treat it independently,
ignoring the occupation of the other levels. The SCF method uses the Fermi function
assuming that the energy of each level depends on the occupation of the other levels. But
this is only approximate. The exact method is to abandon the Fermi function altogether
and use Eq. (3.4.5) instead to calculate the probabilities of the different multi-particle
states.
One well-known example of this is the fact that localized donor or acceptor levels
(which have large charging energies U0 ) in semiconductors at equilibrium are occupied
according to a modified Fermi function (ν is the level degeneracy)
1
f = (3.4.7)
1 + (1/ν) exp [(ε − µ)/kB T ]
rather than the standard Fermi function (cf. Eq. (1.1.1)). We can easily derive this
relation for two spin-degenerate levels (ν = 2) if we assume that the charging energy
U0 is so large that the 11 state has zero probability. We can then write for the remaining
states
α Nα Eα Pα
00 0 0 1/Z
01 1 ε (1/Z ) exp[(µ − ε)/kB T ]
10 1 ε (1/Z ) exp[(µ − ε)/kB T ]
so that Z = 1 + 2 exp[(µ − ε)/kB T ] and the average number of electrons is given by

2 exp[(µ − ε)/kB T ]
N = Nα Pα = P01 + P10 =
α 1 + 2 exp[(µ − ε)/kB T ]
1
=
1 + (1/2) exp[(ε − µ)/kB T ]
in agreement with Eq. (3.4.7). This result, known to every device engineer, could thus
be viewed as a special case of the general result in Eq. (3.4.5).
Equation (3.4.5), however, can only be used to treat equilibrium problems. Our
primary interest is in calculating the current under non-equilibrium conditions and that
is one reason we have emphasized the master equation approach based on Eq. (3.4.2).
For equilibrium problems, it gives the same answer. However, it also helps to bring out
an important conceptual point. One often hears concerns that the law of equilibrium
is a statistical one that can only be applied to large systems. But it is apparent from
the master equation approach that the law of equilibrium (Eq. (3.4.5)) is not a property
of the system. It is a property of the contacts or the “reservoir.” The only assumptions
we have made relate to the energy distribution of the electrons that come in from the
contacts. As long as these “reservoirs” are simple, it does not matter how complicated
or how small the “system” is.
Current calculation: Getting back to non-equilibrium problems, once we have solved

the master equation for the individual probabilities, the source current can be obtained
from

I1 = −q (±) R1 (α → β) Pα
β
+ if b has one more electron than α
− if b has one less electron than α
where R1 represents the part of the total transition rate R associated with the source
contact. In our present problem this reduces to evaluating the expression
- (2γ1 f P00 − γ1 (1 − f )(P01 + P10 )

I1 = ( − q/h) 1 1
+ γ1 f 1 (P01 + P10 ) − 2γ1 (1 − f 1 )P11 ) (3.4.8)
Figure 3.4.3 shows the current–drain voltage (I–VD ) characteristics calculated from the
approach just described. The result is compared with a calculation based on the restricted
SCF method described in Section 1.4. The SCF current–voltage characteristics look
different from Fig. 1.4.6a because the self-consistent potential U0 (N − N0 ) has N0 = 1
rather than zero and we have now included two spins. The two approaches agree well for
U0 = 0.025 eV, but differ appreciably for U0 = 0.25 eV, showing evidence for Coulomb
blockade or single-electron charging (see Exercise E.3.6).
The multi-electron master equation provides a suitable framework for the analysis
of current flow in the Coulomb blockade regime where the single-electron charg-
ing energy U0 is well in excess of the level broadening γ1,2 and/or the thermal
energy kB T . We cannot use this method more generally for two reasons. Firstly, the
(a) U 0 = 0.25 eV (b) U 0 = 0.025 eV

−6
x 10 −6
1.4
1.4 x 10
1.2 1.2
1 1
Current (A)
Current (A)
0.8 0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0 0.2 0.4 0.6 0.8 1
0
0 0.2 0.4 0.6 0.8 1 Drain Voltage, VD (volts)
Drain Voltage, VD (volts)
Fig. 3.4.3 Current vs. drain voltage VD calculated assuming VG = 0 with µ = 0, ε = 0.2 eV, kB T =
0.025 eV, γ 1 = γ 2 = 0.005 eV, UL = −qVD / 2. The two approaches (the SCF and the
multi-electron master equation) agree well for U0 = 0.1 eV, but differ appreciably for U0 = 0.25 eV,
showing evidence for Coulomb blockade or single-electron charging.
size of the problem increases exponentially and becomes prohibitive. Secondly, it is

not clear how to incorporate broadening into this picture and apply it to the trans-
port regime where the broadening is comparable to the other energy scales. And
so it remains a major challenge to provide a proper theoretical description of the
intermediate transport regime U0 ∼ γ1,2 , kB T : the regime where electronic motion
is “strongly correlated” making a two-electron probability like P(11) very different
from the product of one-electron probabilities like P(01)P(10). A lot of work has
gone into trying to discover a suitable SCF within the one-electron picture that will
capture the essential physics of correlation. For example, the self-consistent potential
USCF = U0 N we have used is the same for all energy levels or orbitals. One could use
an “unrestricted” self-consistent field that is orbital-dependent such that the potential
felt by level j excludes any self-interaction due to the number of electrons n j in that
level:
USCF ( j) = U0 (N − n j ) (3.4.9)
Such approaches can lead to better agreement with the results from the multi-electron
picture but must be carefully evaluated, especially for non-equilibrium problems.
EXERCISES
E.3.1. Use the SCF method (only the Hartree term) to calculate the energy of the 1s level
in a helium atom. (a) Plot the nuclear potential UN (r) and the self-consistent electronic
79 Exercises
potential USCF (r) (cf. Fig. 3.1.4a). (b) Plot the wavefunction for the 1s level in helium
and compare with that for the 1s level in hydrogen (cf. Fig. 3.1.4b).
E.3.2. Use the SCF method (only the Hartree term) to calculate the energies of the 3s
and 3p levels in a silicon atom. Plot the wavefunction for the 1s and 3p levels in silicon
and compare with that for the 1s level in hydrogen (cf. Fig. 3.1.4b).
E.3.3. Plot the approximate binding energy for a hydrogen molecule as a function of
the hydrogen–hydrogen bond length, making use of Eqs. (3.3.6) and (3.3.9a, b) and
compare with Fig. 3.3.4.
E.3.4: In Section 1.2 we obtained the following expression for the current through a
single level
q γ1 γ2
I = [ f 1 (ε) − f 2 (ε)]
h γ1 + γ2
and for the average number of electrons
γ1 f 1 + γ2 f 2
N=
γ1 + γ2
by writing a set of rate equations for a single one-electron energy level (without spin
degeneracy). In the multi-electron picture we have two levels “0” and “1” corresponding
to the one-electron level being empty or full respectively. Write down the appropriate
rate equations in this picture and re-derive the expressions for “N” and “I”.
E.3.5: Consider a channel with two spin-degenerate levels assuming the following
parameters: µ = 0, ε = 0.2 eV, kB T = 0.025 eV, γ1 = γ2 = 0.005 eV.
(a) Calculate the number of electrons vs. gate voltage VG , with VD = 0 and UL =
−q VG , using (1) the multi-electron master equation and (2) a restricted SCF
potential given by USCF = U0 (N − N0 ) with N0 = 1. Use two different values
of U0 = 0.025 eV, 0.25 eV.
(b) Calculate the current vs. drain voltage VD assuming VG = 0 with UL = −q VD /2,
using (1) the multi-electron master equation and (2) the restricted SCF potential
given in (a). Use two different values of U0 = 0.025 eV, 0.25 eV and compare with
Fig. 3.4.3.
(c) Repeat (a) and (b) with an unrestricted SCF potential (Eq. (3.4.9)) that excludes
the self-interaction:
Uscf (↑) = U0 (N − n ↑ ) = U0 (n ↓ ) = U0 (n ↓ − 0.5)

Uscf (↓) = U0 (n ↑ − 0.5)
Note: The result may be different depending on whether the initial guess is symmetric,
Uscf (↑) = Uscf (↓) or not, Uscf (↑) = Uscf (↓).
E.3.6: In Fig. 3.4.3a (U0 = 0.25 eV) the multi-electron approach yields two current
plateaus: a lower one with ε2 + UL > µ1 > ε1 + UL such that f 1 1, f 1 0 and
an upper one with µ1 > ε2 + UL > ε1 + UL , such that f 1 1, f 1 1. In either case
f 2 0, f 2 0. Show from Eqs. (3.4.3) and (3.4.8) that the current at these plateaus
is given by
2γ1 γ2 2γ1 γ2
and
2γ1 + γ2 γ1 + γ2
respectively.

CH 2

Uploaded by

Copyright:

Available Formats

CH 2

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

CH 2

Uploaded by

Copyright:

Available Formats

2 Schrödinger equation

2.1 Hydrogen atom

2πr = n(h/mv) (2.1.3)

rn = (n 2 /Z )a0(Bohr radius) (2.1.4)

(a) Ground state (b) Excited states

x=0 x=L x=0 x=L

u = u 0 exp(ikx) exp(−iωt) ⇒ ω2 = v2 k 2 (2.1.11)

u = u 0 sin(kx) exp(−iωt) ⇒ k = nπ/L ⇒ ω = nπ v/L (2.1.12)

Fig. 2.1.4 Standing waves in a ring.

Electron density and probability current density: An electron with a wavefunction

2.2 Method of finite differences

{ψ1 ψ2 ··· · · ·} = {(x1 ) (x2 ) ··· · · ·}

U (x) (x) → U (xn )(xn )

This allows us to write (note: t0 ≡ h- 2 /2ma 2 and Un ≡ U (xn ))

= [(Un + 2t0 ) δn,m − t0 δn,m+1 − t0 δn,m−1 ]ψm (2.2.1)

The elements of the Hamiltonian matrix are given by

Hn,m = [Un + 2t0 ] δn,m − t0 δn,m+1 − t0 δn,m−1 (2.2.3)

[H ] {α} = E α {α} (2.2.5)

2.3.1 Particle in a box

φα (x) ∼ sin(kα x) where kα = απ/L , α = 1, 2, . . .

and their energies are given by E α = h- 2 kα2 /2m.

Normalization: Figure 2.3.2b shows the eigenfunction squared corresponding to the

Numerically, a normalized eigenvector satisfies the condition

So when we compare numerical results with analytical results we should expect

|φα (xn )|2 = |φα (x)|2 a (2.3.2)

−t0 ψ0 + (2t0 + U1 )ψ1 − t0 ψ2 with (2t0 + U1 )ψ1 − t0 ψ2

−t0 ψ N −1 + (2t0 + U N )ψ N − t0 ψ N +1 with −t0 ψ N −1 + (2t0 + U N )ψ N

Mathematically, periodic boundary conditions are implemented by modifying the

φα (x) ∼ sin(kα x) and cos(kα x)

where kα = α 2π/L , α = 1, 2, . . . instead of

As we discussed earlier, instead of writing the eigenstates as

cos(kα x) and sin(kα x)

we could just as well write them as

eikα x = cos(kα x) + i sin(kα x) and e−ikα x = cos(kα x) − i sin(kα x)

sin(kα x) and sin([kα + (2π/a)]x)

represent the same state because at any lattice point xn = na,

sin(kα xn ) = sin([kα + (2π/a)]xn )

2.3.2 Particle in a 3D “box”

Spherically symmetric potential: Some problems may not be separable in Cartesian

(r, θ, φ) = [ f (r )/r ] Ylm (θ,φ) (2.3.7)

where the radial wavefunction f (r ) is obtained by solving the radial Schrödinger

Since, from Eq. (2.3.7)

(r, θ, φ) = [f (r )/r ]Ylm (θ,φ)

and the spherical harmonics are normalized such that

E 1s = −13.56 eV and E 2s = −2.96 eV

assuming γ is (a) purely imaginary (= iβ) and (b) purely real.

3.1 The self-consistent field (SCF) procedure

r ) = −2q 2 /4π ε0r

r ) = −Zq 2 /4π ε0r we will obtain energy

E n = −(Z 2 /n 2 ) E 0 = −54.4 eV/n 2 (Z = 2)

He+ + hν → He++ + e− (3.1.1b)

Step 1. Guess electronic potential USCF (

(a) Nuclear charge Electronic charge

First term in Second term in

If we assume the charge distribution to be isotropic, we can write

For Step 4 we can use straightforward electrostatics to show that

We obtain the total potential by adding the two components.

3.2 Relation to the multi-electron picture

Multi-electron Schrödinger equation: It is important to recognize that the SCF

{ψ1 ψ2 ··· · · ·} = {(x1 ) (x2 ) ··· · · ·}

U (x) (x) → U (xn )(xn )

(r, θ, φ) = [ f (r )/r ] Ylm (θ,φ) (2.3.7)

(r, θ, φ) = [f (r )/r ]Ylm (θ,φ)

E(Na+ ) + E(Cl− ) − E(Na+ Cl− ) = U0

E bin = E(Na) + E(Cl) − E(Na+ Cl− ) = 8.8 eV − U0 + U0 (3.3.3)

E(2H) = Ue,N + Ue ,N = 2E 0 (3.3.7a)

E(H2 ) = UN,N + Ue,e + Ue,N + Ue,N + Ue ,N + Ue ,N (3.3.7b)

E(H2 ) = UN,N + Ue,e + 2E B0 (3.3.7c)

E bin = E(H2 ) − E(2H) = 2(E B0 − E 0 ) + UN,N + Ue,e (3.3.8)

UN,N = q 2 /4π ε0 R (3.3.9a)

E B = Ue,e + Ue,N + Ue,N = Ue,e + E B0 (3.3.10)

E bin = 2(E B − E 0 ) + UN,N

since EB includes the electron–electron interaction energy Ue,e . However, it is easy to

E bin = 2(E B − E 0 ) + (UN,N − Ue,e )

f 1 ≡ f 0 (ε1 − µ1 ) and f 2 ≡ f 0 (ε1 − µ2 )

f 1 ≡ f 0 (ε2 − µ1 ) and f 2 ≡ f 0 (ε2 − µ2 )

f 1 = f 0 (ε1 + UL − µ1 ), f 2 = f 0 (ε1 + UL − µ2 ) (3.4.4a)

- (2γ1 f P00 − γ1 (1 − f )(P01 + P10 )

USCF ( j) = U0 (N − n j ) (3.4.9)

Uscf (↑) = U0 (N − n ↑ ) = U0 (n ↓ ) = U0 (n ↓ − 0.5)