7_tight_binding.md

import matplotlib
from matplotlib import pyplot

import numpy as np

from common import draw_classic_axes, configure_plotting

configure_plotting()

pi = np.pi

Electrons and phonons in 1D

(based on chapters 9–11 of the book)

!!! summary "Learning goals"

After this lecture you will be able to:

- formulate equations of motion for electrons and phonons in 1D.
- solve these equations to arrive at the corresponding dispersion relations.
- derive the group velocity, effective mass, and density of states from the dispersion relation.

Last lecture:

• Movement of a few atoms (Newton's equations)
• Energies of molecular orbitals (LCAO)

This lecture:

• Infinitely many atoms
• Main idea: use periodicity in space, similar to periodicity in time

Equations of motion

Phonons

If atom n has displacement u_n, then

m \ddot{u}_n = -\kappa (u_n - u_{n-1}) -\kappa (u_n - u_{n+1}),

and if we have a system of size L = Na with periodic boundary conditions, we additionally have u_L = u_0.

Electrons

A similar problem: we consider an LCAO wave function \psi = \sum_n c_n \phi_n. We need to solve the Schrödinger equation

E c_n = E_0 c_n + t c_{n+1} + t c_{n-1}

The periodic boundary conditions mean c_N = c_0.

Key idea

Our task is to find all normal modes and to take L → ∞.

All atoms are similar ⇒ the solution should be similar for all atoms ⇒
For phonons: u_n = e^{i \omega t - i k x} u_0, and for electrons: c_n = e^{i E t/\hbar - i k x} c_0. (Here we wrote the time-dependent solution of the Schrödinger equation to emphasize the similarity between two systems.)

Periodicity requires e^{i k L} = 1 ⇒ k = 2\pi p/L, with p\in \mathbb{Z}.

In that case c_n/c_0 = \exp(2 \pi n p a/L) = \exp(2 \pi n p/N), and changing p→p+N corresponds to exactly the same solution, and we have N different solutions in total.

We can see that the solutions with different k (or different p) are identical by plotting them:

x = np.linspace(-.2, 2.8, 500)
fig, ax = pyplot.subplots()
ax.plot(x, np.cos(pi*(x - .3)), label=r'$k=\pi/a$')
ax.plot(x, np.cos(3*pi*(x - .3)), label=r'$k=3\pi/a$')
ax.plot(x, np.cos(5*pi*(x - .3)), label=r'$k=5\pi/a$')
sites = np.arange(3) + .3
ax.scatter(sites, np.cos(pi*(sites - .3)), c='k', s=64, zorder=5)
ax.set_xlabel('$x$')
ax.set_ylabel('$u_n$')
ax.set_xlim((-.1, 3.2))
ax.set_ylim((-1.3, 1.3))
ax.legend(loc='lower right')
draw_classic_axes(ax)
ax.annotate(s='', xy=(.3, -1.1), xytext=(1.3, -1.1),
            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
ax.text(.3 + .5, -1.25, '$a$', ha='center');

Let's count how many different solutions we expect to find. We have a system with N degrees of freedom (either u_n or c_n), and therefore we need N normal modes (or eigenstates).

Because we proposed N different plane wave solutions, if we find an energy or frequency for each solution, we have fully solved the problem!

Solution

Phonons

First substitute the Ansatz into the equations of motion: -m \omega^2 c_0 e^{i\omega t - ikx} = \kappa c_0 e^{i\omega t}(-2 e^{-ikx} + e^{-ikx+ika}+ e^{-ikx-ika}), the exponents and c_0 drop out (just like we expected) and we get: -m \omega^2 = \kappa (-2 + e^{ika}+ e^{-ika})=\kappa [-2 + 2\cos(ka)], or after a further simplification: \omega = \sqrt{\frac{2\kappa}{m}}\sqrt{1-\cos(ka)} = 2\sqrt{\frac{\kappa}{m}}|\sin(ka/2)| [here we used that 1-\cos(x)=2\sin^2(x/2).]

So we arrive to the dispersion relation

k = np.linspace(-2*pi, 6*pi, 500)
fig, ax = pyplot.subplots()

pyplot.plot(k, np.abs(np.sin(k/2)))

ax.set_ylim(bottom=0, top=1.2)
ax.set_xlabel('$ka$')
ax.set_ylabel(r'$\omega$')
pyplot.xticks(list(2*pi*np.arange(-1, 4)) + [-pi, pi],
              [r'$-2\pi$', '$0$', r'$2\pi$', r'$4\pi$', r'$6\pi$',
               r'$-\pi$', r'$\pi$'])
pyplot.yticks([1], [r'$2\sqrt{\frac{\kappa}{m}}$'])
pyplot.vlines([-pi, pi], 0, 1.1, linestyles='dashed')
pyplot.hlines([1], .1*pi, 1.3*pi, linestyles='dashed')
draw_classic_axes(ax)
ax.annotate(s='', xy=(-pi, -.15), xytext=(pi, -.15),
            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
ax.text(0, -.25, '1st Brillouin zone', ha='center')
ax.set_ylim(bottom=-.7);

Here k_p = 2\pi p / L for 0 ≤ p < N are all the normal modes of the crystal, and therefore we can rederive a better Debye model!

Before we had \sum_p → \frac{L}{2\pi}\int_{-\omega_D/v_s}^{\omega_D/v_s}dk, because we introduced the cutoff to fix the problem with the number of modes.

Now \sum_p → \frac{L}{2\pi}\int_{-\pi/a}^{\pi/a}dk, the integral is over a finite number of modes because there is only a finite number of modes.

Sound velocity: at small k \sin(ka/2)\approx ka/2, and therefore \omega \approx \sqrt{\kappa/m} k a, so we have derived the existence of sound waves.

Electrons

Substitute the Ansatz into the equations of motion:

E c_0 e^{-ikx} = E_0 c_0 e^{-ikx} +t e^{-ikx-ika} + te^{-ikx+ika}, and after canceling the exponents we immediately get E = E_0 + 2t\cos(ka), so we come to the dispersion relation:

pyplot.figure()
k = np.linspace(-pi, pi, 300)
pyplot.plot(k, -np.cos(k))
pyplot.xlabel('$ka$'); pyplot.ylabel('$E$')
pyplot.xticks([-pi, 0, pi], [r'$-\pi$', 0, r'$\pi$'])
pyplot.yticks([-1, 0, 1], ['$E_0+2t$', '$E_0$', '$E_0-2t$']);

Usually electron dispersion has multiple options for E(k), each called a band. The complete dispersion relation is also called a band structure.

The band bottom is at k=0 (when t<0). There the energy

E = E_0 + 2t\cos{ka} \approx E_0 + 2t - t (ka)^2. Comparing this with E=(\hbar k)^2/2m, we see that the dispersion relation is as if the electrons had an effective mass m^*=\hbar^2/2ta^2.

A single band has 2N states (N k-values and 2 spins). If each atom contributes 2 electrons, all these states must be occupied. Applying electric field, therefore cannot create electric current ⇒ we have explained insulators!

Group velocity, effective mass, density of states

(here I only discuss electrons; for phonons everything is the same except for replacing E = \hbar \omega)

Question: what happens if we apply an external electric field to the crystal:

A reminder from Hamiltonian mechanics:

H = \frac{p^2}{2m} + U(r)

\begin{aligned} \frac{dr}{dt} &= \frac{\partial H(p, r)}{\partial p} \equiv v \\ \frac{dp}{dt} &= -\frac{\partial H(p, r)}{\partial r} \equiv F \end{aligned}

The first equation is velocity, the second equation is force.

Connection to quantum mechanics: p = \hbar k.

If electrons form bands, then in each band H = E(k) + \tilde{U}(r),

where \tilde{U}(r) = -e\mathbf{E} only includes slow variations of the electrostatic potential (the rapidly changing atomic potential is responsible for E(k)):

For electrons in an arbitrary band structure: \hbar^{-1}\partial H/\partial k is the group velocity (same as for phonons).

Another important concept is the effective mass defined using an analogy:

\frac{dv}{dt} = m_{eff}^{-1} F.

Substituting, we get

m_{eff} = \hbar^2\left(\frac{d^2 E(k)}{dk^2}\right)^{-1}

pyplot.figure(figsize=(8, 5))
k = np.linspace(-pi, pi, 300)
meff = 1/np.cos(k)
color = list(matplotlib.rcParams['axes.prop_cycle'])[0]['color']
pyplot.plot(k[meff > 0], meff[meff > 0], c=color)
pyplot.plot(k[meff < 0], meff[meff < 0], c=color)
pyplot.ylim(-5, 5)
pyplot.xlabel('$ka$'); pyplot.ylabel('$m_{eff}$')
pyplot.xticks([-pi, 0, pi], [r'$-\pi$', 0, r'$\pi$']);

density of states

DOS is the number of states per unit energy. In 1D we have g(E) = \sum |dk/dE| = \sum |v|^{-1},

The sum goes over all bands existing at a given energy, positive and negative momenta, and spins.

Using

E = E_0 + 2t \cos(ka),

we get

k = \pm\arccos[(E - E_0) / 2t], and |d k / d E| = [4t^2 - (E - E_0)^2]^{-1/2}.

You can get to this result immediately if you remember the derivative of arccosine. Otherwise you need to go the long way: compute dE/dk as a function of k, express k through E as we did above, and take the inverse.

Also a sanity check: when the energy is close to the bottom of the band, E = E_0 + 2t + \delta E we get g(E) \sim E^{-1/2}, as we would expect in 1D.