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 more than one degree 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 hoppings in LCAO

Consider electrons in a 1D atomic chain again

Let's relax the assumption that there is no hopping between atoms that are further than nearest neighbors:

$$\langle \phi_n | H | \phi_{n+2}\rangle \equiv t' ≠ 0.$$

We now have more terms in the Schrödinger equation:

$$E c_n = E_0 c_n + t c_{n-1} + t c_{n+1} + t' c_{n-2} + t' c_{n+2},$$

but the same Ansatz as before should work: $ c_n = c_0 \exp(ikna) $.

Indeed, after substituting we get:

$$E = E_0 + 2t\cos(ka) + 2t'\cos(2ka).$$

(check that!)

Therefore including more hopping terms works out of the box using the same Ansatz.

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 whem to efficiently emit and adsorb photons.

Important concept: density of states (DOS):

$g(\omega)=\rho_{\rm S}\frac{ {\rm d}k}{ {\rm d}\omega}$

The quantity

$\frac{ {\rm d}k}{ {\rm d}\omega}$

can be calculated from the dispersion.

When plotted, the DOS may look something like this:

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 normally

$g(\omega)$

is 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$

? 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 reconsile the two pictures let's plot two unit cells in the 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)

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}...$

$\rightarrow \rho_{\rm S}(k)=$

(no. of

$k$

-values / unit of

$k$

)

$=\frac{L}{\pi}$

, where only the positive part of

$k$

-space is filled with allowed

$k$

-values.

For open boundary conditions

$\rightarrow$

running waves

$\rightarrow k=$

$...\frac{-4\pi}{L},\frac{-2\pi}{L},0,\frac{2\pi}{L},\frac{4\pi}{L}...$

$\rightarrow \rho_{\rm R}(k)=\frac{L}{2\pi}$

, which is lower than for the case of closed boundary conditions, however, in this case the entire

$k$

-space is filled.

Summary

• By using plane waves in real space as an Ansatz we found all normal modes and eigenvectors
• Dispersion relation of a system with period
  $a$
  in real space is periodic with period
  $2\pi/a$
  in
  $k$
  -space
• Computing dispersions explains all the problems we listed before (need for cutoff, lack of scattering with every single atom, existence of insulators).
• Electrons and phonons have a complicated nonlinear relation between momentum and velocity (group velocity), effective mass, and density of states
• In a system with more than one degree of freedom per unit cell we need to consider independent amplitudes for each degree of freedom, and get multiple bands.