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.1-9.3 & 11.1-11.3 of the book)

!!! success "Expected prior knowledge"

Before the start of this lecture, you should be able to:

- Derive Newton's equations of motion of a triatomic chain (previous lecture).
- Write down the dispersion relation of phonons in the Debye model.
- Express complex exponentials through trigonometric functions and vice versa.
- Apply Taylor expansion trigonometric functions.
- Take derivatives inverse trigonometric functions.

!!! summary "Learning goals"

After this lecture you will be able to:

- Formulate the equations of motion of electrons and phonons in 1D.
- Derive the dispersion relation from the equations of motion.
- Derive the group velocity, effective mass, and density of states from the dispersion relation.

Last lecture:

• Vibrational modes of few-atom chains (analyzed using Newton's equations)
• Orbital energies of electrons in few-atom chains (analyzed using LCAOs)

This lecture:

• Phonons and electrons in chains of infinitely many atoms.
• Main idea: use periodicity in space, similar to periodicity in time

To emphasize the similarities and the differences between electrons and phonons, we will deal with both types of particles at once.

Equations of motion

Phonons

In the Debye model, we assumed that the dispersion relation is strictly linear in k. Now is the time to revisit this assumption. To do that, let us consider a 1D homogeneous chain of atoms. We assume that the atoms in the chain interact only with their nearest neighbors through a harmonic potential, like we derived in the previous lecture. In other words, we model the atoms as point masses connected by identical springs.

We denote the displacement of atom n from equilibrium by u_n. Within this convention, Newton's equation of motion for the n-th atom is given by:

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

We use the periodic boundary conditions just like we did in the Sommerfield model. The boundary conditions imply that in a system of size L = Na, we have u_N = u_0.

Electrons

The following figure shows the interatomic potential with atoms placed at x_k = ka with k \in \mathbb{Z}:

x = np.linspace(0.001, 3, 200)

fig, ax = pyplot.subplots(1, 1)
ax.plot(x, 1.2-1/np.abs(np.sin(np.pi * x))**(1/2))
ax.set_ylim(-.7, .5)
ax.set_xlabel("$x$")
ax.set_ylabel("$U(x)$")
ax.set_xticks([-.05, 1, 2])
ax.set_xticklabels(["$0$", "$a$", "$2a$"])

draw_classic_axes(ax)

Similarly to the triatomic system case, we formulate the molecular orbital via the LCAO model: \vert \Psi \rangle = \sum_n \phi_n |n \rangle. We assume only nearest-neighbor hopping -t and an on-site energy E_0. The coupled Schrödinger equation of the |n \rangle orbital is: E \phi_n = E_0 \phi_n - t \phi_{n+1} - t \phi_{n-1}.

Again, the periodic boundary conditions imply \phi_N = \phi_0.

We now have the equations of motion of both phonons and electrons. All that remains is to solve them.

Key idea for solving these equations

In order to solve the equations of motion, we need to come up with a reasonable guess. If we take a look at the equations of motion, we see that they are the same for all atoms. To be specific, the structure of the equations is the same no matter what value of n we choose. Since these equations define the solutions, we reason that the solutions should also be independent of the choice of n. As a result, we assume a plane wave solution, also called a plane wave ansatz, with the same amplitude for each atom. In the case of phonons, we obtain

u_n = Ae^{i \omega t - i k x_n},

and the ansatz for electrons is given similarly by

\phi_n = Be^{i E t/\hbar - i k x_n}.

We wrote x_n=na and we wrote the time-dependent solution of the Schrödinger equation to emphasize the similarity between the two systems.

We already know that the periodic boundary conditions only allow plane waves with k being a multiple of 2\pi/L. In the case of the electron system, periodic boundary conditions give \phi_0 = \phi_N, which results in

1 = e^{ik0} = e^{ikNa}.

The above equation defines the allowed values of k:

k = \frac{2 \pi p}{Na}, \quad \text{with} p \in \mathbb{Z}.

We use the quantized values of k in our plane wave ansatz: e^{ikx_n} = e^{i p \frac{2\pi}{Na} n a} = e^{i \frac{2 \pi n p}{N}}. We notice something interesting if we investigate the case of p→p+N. In this case, the plane wave ansatz becomes e^{i \frac{2\pi n(p+N)}{N}} = e^{i \frac{2\pi np}{N} + i2\pi n} = e^{i \frac{2\pi np}{N}}, which is exactly the same solution. Counting the number of inequivalent plane waves, we find exactly N different solutions in total. All that is left is to find the energy of each solution!

The reason why solutions with different values of k are identical is aliasing: because the plane wave is only defined at discrete positions, acquiring a phase factor of 2π between two atoms is equivalent to nothing happening:

x = np.linspace(-.2, 2.8, 500)
fig, ax = pyplot.subplots()
ax.plot(x, np.cos(pi*(x)), label=r'$k=\pi/a$')
ax.plot(x, np.cos(3*pi*(x)), label=r'$k=3\pi/a$')
ax.plot(x, np.cos(5*pi*(x)), label=r'$k=5\pi/a$')
sites = np.arange(3)
ax.scatter(sites, np.cos(pi*(sites)), 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=(0, -1.1), xytext=(1, -1.1),
            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
ax.text(.5, -1.25, '$a$', ha='center');

??? Question "How many different solutions did we expect to find?" We have a system with N degrees of freedom (either u_n or \phi_n), and therefore we expect N normal modes (or eigenstates).

Solving the equations of motion

Phonons

We substitute the plane wave ansatz into the equations of motion:

-m \omega^2 A e^{i\omega t - ikx_n} = \kappa A e^{i\omega t}(-2 e^{-ikx_n} + e^{-ikx_n+ika}+ e^{-ikx_n-ika}).

Searching for solutions with A \neq 0 we obtain

-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)|,

where we substituted 1-\cos(x) = 2\sin^2(x/2).

We arrive at the phonon dispersion relation shown below.

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=-.3);