8_many_atoms.md

import matplotlib
from matplotlib import pyplot

import numpy as np

from common import draw_classic_axes, configure_plotting

configure_plotting()

pi = np.pi

More complex systems

!!! summary "Learning goals"

After this lecture you will be able to:

- formulate equations of motion for electrons or phonons in 1D, for systems with multiple degrees of freedom per unit cell.
- solve these equations to arrive at the dispersion relation.
- derive the group velocity, effective mass, and density of states.
- explain what happens with the band structure when the periodicity of the lattice is increased or reduced.

More degrees of freedom per unit cell:

Before we used that all atoms are similar, but this doesn't work anymore. We need to generalize our anzatz. We label all degrees of freedom in each unit cell (a repeated element of the system).

For two atoms in a unit cell we have displacements u_{1,n} and u_{2,n}. Here the first index is the atom number within the unit cell, second is the unit cell number. The atom with mass M is the atom 1, and the atom with mass m is the atom 2.

Having the degrees of freedom let's write down the equations of motion:

\begin{aligned} M\ddot{u}_{1,n}=\kappa(u_{2,n}-2u_{1,n}+u_{2,n-1})\\ m\ddot{u}_{2,n}=\kappa(u_{1, n} - 2u_{2,n}+u_{1,n+1}), \end{aligned}

The new Ansatz is conceptually the same as before: all unit cells should be the same up to a phase factor:

\begin{pmatrix} u_{1,n}\\ u_{2,n} \end{pmatrix} = e^{i\omega t - ik na} \begin{pmatrix} u_{1}\\ u_{2} \end{pmatrix}.

Substituting this ansatz into the equations of motion we end up with an eigenvalue problem: \omega^2 \begin{pmatrix} M & 0 \\ 0 & m \end{pmatrix} \begin{pmatrix} u_{1} \\ u_{2} \end{pmatrix} = \kappa \begin{pmatrix} 2 & -1 - e^{ika} \\ -1-e^{-ika} & 2 \end{pmatrix} \begin{pmatrix} u_{1}\\ u_{2} \end{pmatrix}, with eigenfrequencies \omega^2=\frac{\kappa(M+m)}{Mm}\pm \kappa\left\{\left(\frac{M+m}{Mm}\right)^2-\frac{4}{Mm}{\rm sin}^2\left(\frac{1}{2}ka\right)\right\}^{\frac{1}{2}}

Looking at the eigenvectors we see that for every k there are now two values of \omega: one corresponding to in-phase motion (–) and anti-phase (+).

def dispersion_2m(k, kappa=1, M=1.4, m=1, acoustic=True):
    Mm = M*m
    m_harm = (M + m) / Mm
    root = kappa * np.sqrt(m_harm**2 - 4*np.sin(k/2)**2 / Mm)
    if acoustic:
        root *= -1
    return np.sqrt(kappa*m_harm + root)

# TODO: Add panels with eigenvectors
k = np.linspace(-2*pi, 6*pi, 300)
fig, ax = pyplot.subplots()
ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
ax.plot(k, dispersion_2m(k), label='acoustic')
ax.set_xlabel('$ka$')
ax.set_ylabel(r'$\omega$')
pyplot.xticks([-pi, 0, pi], [r'$-\pi$', '$0$', r'$\pi$'])
pyplot.yticks([], [])
pyplot.vlines([-pi, pi], 0, 2.2, linestyles='dashed')
pyplot.legend()
pyplot.xlim(-1.75*pi, 3.5*pi)
pyplot.ylim(bottom=0)
draw_classic_axes(ax, xlabeloffset=.2)

The lower branch is called acoustic because its linear dispersion near \omega=0 matches the behavior of the regular sound waves. The upper branch is the optical branch because the high phonon frequency allows them to efficiently emit and adsorb photons.

The density of states (DOS) is given by g(\omega)=\rho_{\rm S}\frac{ {\rm d}k}{ {\rm d}\omega}.

We can derive \frac{ {\rm d}k}{ {\rm d}\omega} from the dispersion relation \omega(k), as graphically shown here:

matplotlib.rcParams['font.size'] = 24
k = np.linspace(-.25*pi, 1.5*pi, 300)
k_dos = np.linspace(0, pi, 20)
fig, (ax, ax2) = pyplot.subplots(ncols=2, sharey=True, figsize=(10, 5))
ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
ax.plot(k, dispersion_2m(k), label='acoustic')
ax.vlines(k_dos, 0, dispersion_2m(k_dos, acoustic=False),
          colors=(0.5, 0.5, 0.5, 0.5))
ax.hlines(
    np.hstack((dispersion_2m(k_dos, acoustic=False), dispersion_2m(k_dos))),
    np.hstack((k_dos, k_dos)),
    1.8*pi,
    colors=(0.5, 0.5, 0.5, 0.5)
)
ax.set_xlabel('$ka$')
ax.set_ylabel(r'$\omega$')
ax.set_xticks([0, pi])
ax.set_xticklabels(['$0$', r'$\pi$'])
ax.set_yticks([], [])
ax.set_xlim(-pi/4, 2*pi)
ax.set_ylim((0, dispersion_2m(0, acoustic=False) + .2))
draw_classic_axes(ax, xlabeloffset=.2)

k = np.linspace(0, pi, 1000)
omegas = np.hstack((
    dispersion_2m(k, acoustic=False), dispersion_2m(k)
))
ax2.hist(omegas, orientation='horizontal', bins=75)
ax2.set_xlabel(r'$g(\omega)$')
ax2.set_ylabel(r'$\omega$')

# Truncate the singularity in the DOS
max_x = ax2.get_xlim()[1]
ax2.set_xlim((0, max_x/2))
draw_classic_axes(ax2, xlabeloffset=.1)
matplotlib.rcParams['font.size'] = 16

Note that g(\omega) is generally plotted along the vertical axis and \omega along the horizontal axis – the plot above is just to demonstrate the relation between the dispersion and the DOS. The singularities in g(\omega) at the bottom and top of each branch are called van Hove singularities.

Consistency check with 1 atom per cell

But what if we take M\rightarrow m? Then we should reproduce the previous result with only one band!

The reason why we get two bands is because we have a\rightarrow 2a compared to 1 atom per unit cell!

For M\rightarrow m, the dispersion diagram should come back to the diagram obtained for the first example (i.e. with identical masses).

To reconcile the two pictures, let's plot two unit cells in reciprocal space. (Definition: each of these unit cells is called a Brillouin zone, a concept that will come up later.)

k = np.linspace(0, 2*pi, 300)
k_dos = np.linspace(0, pi, 20)
fig, ax = pyplot.subplots()
ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
ax.plot(k, dispersion_2m(k), label='acoustic')
omega_max = dispersion_2m(0, acoustic=False)
ax.plot(k, omega_max * np.sin(k/4), label='equal masses')
ax.set_xlabel('$ka$')
ax.set_ylabel(r'$\omega$')
ax.set_xticks([0, pi, 2*pi])
ax.set_xticklabels(['$0$', r'$\pi/2a$', r'$\pi/a$'])
ax.set_yticks([], [])
ax.set_xlim(-pi/8, 2*pi+.4)
ax.set_ylim((0, dispersion_2m(0, acoustic=False) + .2))
ax.legend(loc='lower right')
pyplot.vlines([pi, 2*pi], 0, 2.2, linestyles='dashed')
draw_classic_axes(ax, xlabeloffset=.2)

We see that doubling the lattice constant "folds" the band structure on itself, doubling all the bands.

Total number of states

Which values of k are allowed for a finite 1D crystal of length L?

Assume closed boundary conditions \rightarrow standing waves \rightarrow k=\frac{\pi}{L},\frac{2\pi}{L},\frac{3\pi}{L},\frac{4\pi}{L}...