Many atoms per unit cell

    import matplotlib
from matplotlib import pyplot

import numpy as np

from common import draw_classic_axes, configure_plotting

configure_plotting()

pi = np.pi
  

(based on chapters 10-11 of the book)

Expected prior knowledge

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

• Write down equations of motion and the LCAO Hamiltonian (similar to the previous lectures)
• Solve an eigenvalue problem

Learning goals

After this lecture you will be able to:

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

More degrees of freedom per unit cell:

In the last lecture, we modeled phonons through a 1D homogeneous chain of atoms. The model gave us insight into phonons and justified the crude approximations of the Debye model. Despite the model's usefulness, we are not always able to model problems as a homogeneous chain. Therefore, let us tackle a slightly more general problem. Consider a chain of atoms with two masses, \(m_1\) and \(m_2\). In our chain, an atom with mass \(m_1\) is always followed by an atom of mass \(m_2\) and vice versa. Similar to before, the atoms interact via a harmonic potential with spring constant \(\kappa\).

In the last lecture, we solved the homogenous chain problem by identifying its high level of symmetry. However, the problem here is not as symmetric. We need to re-think and generalize our ansatz.

First we identify a pattern that does repeat: an atom with mass \(m_1\) followed by one with mass \(m_2\). This is called a unit cell (the smallest repeated element of the system)

We now label all degrees of freedom in a unit cell. The two atoms in a unit cell have displacements \(u_{1,n}\) and \(u_{2,n}\), where the first index specifies the atom number within the unit cell and the second the unit cell number. In our choice of unit cell, atom 1 has mass \(m_1\) and atom 2 has mass \(m_2\).

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

\[\begin{align*} m_1\ddot{u}_{1,n}&=\kappa(u_{2,n}-2u_{1,n}+u_{2,n-1})\\ m_2\ddot{u}_{2,n}&=\kappa(u_{1, n} - 2u_{2,n}+u_{1,n+1}). \end{align*}\]

The new ansatz is conceptually the same as before: all unit cells should behave 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} A_{1}\\ A_{2} \end{pmatrix}. \]

Substituting this ansatz into the equations of motion (and assuming that the solution is nontrivial) we end up with an eigenvalue problem:

\[ \omega^2 \begin{pmatrix} m_1 & 0 \\ 0 & m_2 \end{pmatrix} \begin{pmatrix} A_{1} \\ A_{2} \end{pmatrix} = \kappa \begin{pmatrix} 2 & -1 - e^{ika} \\ -1-e^{-ika} & 2 \end{pmatrix} \begin{pmatrix} A_{1}\\ A_{2} \end{pmatrix}, \]

with eigenfrequencies:

\[ \omega^2=\frac{\kappa(m_1+m_2)}{m_1m_2}\pm \kappa\left\{\left(\frac{m_1+m_2}{m_1m_2}\right)^2-\frac{4}{m_1m_2}\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 positive 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)
ax.annotate('', xy=(-pi, -.3), xytext=(pi, -.3),
            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
ax.text(0, -.55, '1st Brillouin zone', ha='center');
#draw_classic_axes(ax, xlabeloffset=.2)

The figure above shows a plot of the eigenfrequencies as a function of \(ka\). Just like last time, the plot is periodic in \(k\), with \(k\)-values which differ by a multiple of \(2 \pi / a\) corresponding to the same solution. Therefore we only look at the \(k\)-values between \(k - \pi / a\) and \(k + \pi / a\). This range (in between the dashed lines) is called the first Brillouin zone, and functions as a unit cell in reciprocal space. The Brillouin zone will be explained in more detail in the lecture about X-ray diffraction.

Unlike the simple atomic chain, the dispersion relation now has two branches (or bands). The reason an additional branch appears in the solution is due to the 2 degrees of freedom per unit cell. If we had started with 3 different atoms, the eigenvalue problem would contain 3 equations, so that there would be three eigenfrequencies per \(k\) and the dispersion relation would have 3 branches.

The lower branch is called acoustic because its linear dispersion near \(\omega=0\) matches the behavior of regular sound waves. The upper branch is the optical branch because it can cross with the (linear) dispersion relation of photons, allowing these phonons to efficiently emit and absorb photons.

Like before, the phonon group velocity is \(v=\textrm{d}\omega/\textrm{d}k\), and the density of states is \(g(\omega)=\textrm{d}N/\textrm{d}\omega = \frac{L}{2π} ∑| \textrm{d}k/\textrm{d} \omega|\). The sum goes over all states at a given energy. In this case, the sum ensures we include the contribution to the DOS from both the positive and negative momenta. Since the energy of this system is symmetric with respect to momentum reversal, the sum only introduces a factor of 2.

An intuitive way to visualize the density of states \(g(\omega)\) is to consider it a histogram of the samples drawn from the dispersion relation \(\omega(k)\):

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_yticklabels([])
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))
ax2.set_xticks([])
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 right 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

To check if our result is consistent with the previous lecture, we examine what happens when we take \(m_1\rightarrow m_2\). Physically, the system is then exactly the same as the 1D monatomic chain, so we would want our solutions to also be the same. But at first glance the two solutions seem very different: the solution with one atom has only one band, but the solution with two atoms has two bands.

To reconcile the two pictures, let's plot two unit cells in reciprocal space.

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('$k$')
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_yticklabels([])
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 must be careful when considering the lattice constant \(a\). In the 1D monatomic chain, \(a\) was the distance between two neighbouring atoms, but in our diatomic chain \(a\) is the size of the unit cell, which is twice as big. Therefore we take the size of our unit cell in the diatomic case to be \(2a\), so the physical systems are the same.

Looking at the graph we see that doubling the lattice constant "folds" the band structure on itself. There are then \(2\) possible values of \(\omega\) for each \(k\), creating the \(2\) branches of the non-periodic system. We can look at the graph in two distinct ways: We can consider the green line (monatomic chain solution). Its Brillouin zone extends from \(0\) to \(\pi / a\). Alternatively, we consider the orange and blue lines together. In this case, the Brillouin zone is smaller ranging from \(0\) to \(\pi / 2a\) (because of the larger lattice constant).

Despite the two band structures look different, the density of states only changes very little: when \(m_1 ≈ m_2\) only the states near the band touching point slightly move in frequency.

Summary

• By using plane waves in real space as an ansatz, we found all normal modes of an atom chain with two different atoms. (Just like in the case of 1 degree of freedom per unit cell).
• The density of states can be derived graphically from the dispersion relation.
• The dispersion relation of a system with period \(a\) in real space is periodic with period \(2\pi/a\) in \(k\)-space.
• 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 we get multiple bands.
• Systems with different band structures can have the same density of states.

Exercises

Warm-up questions

1. Verify that the expression for \(\omega^2\) found in this lecture is always positive. Why is this important?
2. When calculating the DOS, we only look at the first Brillouin zone. Why?
3. Consider a chain of identical atoms with onsite energy \(\varepsilon_0\) connected by hoppings \(-t\). Write down the LCAO wavefunction. Write down the Schrodinger equation for the probability amplitudes. Substitute the Ansatz and solve for the dispersion.
4. Suppose we let the hoppings of the monatomic chain alternate between two slightly different values (see the sketch in the figure of exercise 2), thereby increasing the size of the unit cell. At what value of \(k\) do you expect band gaps to open in the dispersion derived in warmup question 3?

Exercise 1: Analyzing the diatomic vibrating chain

In this lecture, we derived the eigenfreqencies of a diatomic vibrating chain with 2 different, alternating masses. We found that a gap opens in the vibrational spectrum. Here we analyze the group velocity, density of states and band gap of the diatomic chain.

Defining the reduced mass \(\mu = \tfrac{m_1m_2}{m_1+m_2}\), the eigenfrequencies can be written as

\[ \omega^2=\frac{\kappa}{\mu}\left(1 \pm \sqrt{1-\frac{4\mu^2}{m_1m_2}\sin^2\left(\frac{1}{2}ka\right)}\right), \]

where the plus sign corresponds to the optical branch and the minus sign to the acoustic branch.

1. Work out the expression for \(\omega^2\) assuming \(m_1 = m_2\). Compare this to the solution for the monatomic chain.
2. Now take \(m_1\) and \(m_2\) to be slightly different. Sketch the dispersion. Indicate the frequency at which the band gap opens due to the different masses

3. Find the magnitude of the group velocity near \(k=0\) for the acoustic branch.

   Hint

   Make use of a Taylor expansion.

4. Show that the group velocity at \(k=0\) for the optical branch is zero.

5. Derive an expression of the density of states \(g(\omega)\) for the acoustic branch and small \(k\). Make use of your expression of the group velocity. Compare this expression with the 1D Debye density of states derived in exercise 1 of the Debye lecture.
6. Show that the size of the band gap is \(\Delta\omega \approx \sqrt{\frac{\kappa}{\mu}}\frac{|m_1-m_2|}{m_1+m_2}\) when \(m_1\approx m_2\).

Exercise 2: Analyzing the LCAO chain with alternating hoppings

In the previous lecture, we derived the electronic band structure of a 1D, equally-spaced atomic chain, with the equal spacing causing identical nearest-neighbour hoppings throughout the chain. It turns out however that, depending on the number of electrons in the system, such chains are unstable: The reason is that distorting the equal spacing can lower the energy of the system. This structural change is known as the Peierls transition.

The spacing of the distorted chain will alternate between two different distances. This causes the hopping energy to alternate between \(t_1\) and \(t_2\). In this exercise we analyze the band structure of this LCAO chain with alternating hoppings. We set the onsite energies for all atoms to \(\varepsilon_0\). The situation is depicted in the figure below.

Due to the alternating hoppings, we have two atoms in each unit cell. We call the associated atomic orbitals in the n-th unit cell \(|n,1\rangle\) and \(|n,2 \rangle\), as indicated in the figure. We construct the electronic states of the entire chain using linear combinations of the atomic orbitals

\[ \left|\Psi \right\rangle = \sum_n \left(\phi_n \left| n,1 \right\rangle + \psi_n \left| n,2 \right\rangle\right) \]

As usual, we assume the atomic orbitals to be orthogonal to each other.

1. Copy the figure and indicate the unit cell. Then, using the Schrödinger equation, derive the equations of motion describing the electron eigenstates.

   Hint

   To this end, find expressions for \(E \left< n,1 \vert \Psi \right> = \left< n,1 \right| H \left|\Psi \right>\) and \(E \left< n,2 \vert \Psi \right> = \left< n,2 \right| H \left|\Psi \right>\).

2. Using the Ansatz \(\phi_n = \phi_0 e^{ikna}\) and \(\psi_n = \psi_0 e^{ikna}\), show that the Schrödinger equation can be written as:

   \[ \begin{pmatrix} \varepsilon_0 & t_1 + t_2 e^{-i k a} \\ t_1 + t_2 e^{i k a} & \varepsilon_0 \end{pmatrix} \begin{pmatrix} \phi_0 \\ \psi_0 \end{pmatrix} = E \begin{pmatrix} \phi_0 \\ \psi_0 \end{pmatrix}. \]

3. Derive the dispersion relation from this eigenvalue problem and sketch it. Indicate the size of the band gap. Does the dispersion reduce to that of the 1D, equally-spaced atomic chain if \(t_1 = t_2\)? Does it resemble the band structure shown on the Wikipedia page?

4. Sketch the density of states by graphically constructing it from your sketch of the dispersion (no calculations!)

5. Calculate the group velocity \(v(k)\) for both bands. From the group velocity, calculate the density of states \(g(E)\) of the entire band structure and make a plot of it. Check if the plot matches your sketch of subquestion 5.
6. What are the two possible eigenvectors \([\phi_0 \quad \psi_0]\) at \(k=0\)? Use this to sketch the two possible wavefunctions \(|\Psi_\pm(k=0)\rangle\) as a function of the coordinate \(x\) along the chain (assume some peak-shaped atomic orbital such as a Gaussian). Discuss how to understand which wave function has the lowest energy.
7. Repeat the previous subquestion for \(k=\pi/a\). Use your sketch to argue why the two wavefunctions \(|\Psi(k=\pi/a)\rangle\) have different energy when \(t_1 \neq t_2\), and why they are degenerate when the \(t_1=t_2\).
8. Consider your sketch of the density of states. Where is the Fermi energy? Argue why a chain of equidistant atoms can lower its energy by displacing the atoms such that the hoppings become unequal (this is the Peierls transition). Is the same true if each atom contributes 2 electrons?

Exercise 3: Atomic chain with 3 different spring constants

Suppose we have a vibrating 1D atomic chain with 3 different spring constants alternating like \(\kappa_ 1\), \(\kappa_2\), \(\kappa_3\), \(\kappa_1\), etc. All the atoms in the chain have an equal mass \(m\).

1. Make a sketch of this chain and indicate the length of the unit cell \(a\) in this sketch.
2. Derive the equations of motion for this chain.

3. By filling in the trial solutions into the equations of motion (which should be similar to Ansatz used in the lecture), show that the eigenvalue problem is

   \[ \omega^2 \mathbf{u} = \frac{1}{m} \begin{pmatrix} \kappa_1 + \kappa_3 & -\kappa_1 & -\kappa_3 e^{i k a} \\ -\kappa_1 & \kappa_1+\kappa_2 & -\kappa_2 \\ -\kappa_3 e^{-i k a} & -\kappa_2 & \kappa_2 + \kappa_3 \end{pmatrix} \mathbf{u} \]

4. In general, the eigenvalue problem above cannot be solved analytically. It can however be solved in specific cases. Find the eigenvalues \(\omega^2\) when \(ka = 0\) and \(\kappa_1 = \kappa_2\) .

   Hint

   To solve the eigenvalue problem quickly, make use of the symmetry of the problem to guess the form of the eigenvectors (make a drawing!) as we did in exercise 1 of Bonds and Spectra. Alternatively, and more formally, you can use that the mass-spring matrix commutes with the matrix

   \[ X = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{pmatrix}. \]

   Commuting matrices have at least one identical eigenvector, and their eigenvectors that belong to a degenerate eigenvalue span the same eigenspace.

5. Also find the eigenvalues \(\omega^2\) when \(ka = \pi/a\) and \(\kappa_1 = \kappa_2\), again making use of the symmetry of the problem.

6. What will happen to the periodicity of the band structure if \(\kappa_1 = \kappa_2 = \kappa_3\)? Sketch the band structure for this case. Now suppose \(\kappa_3\) becomes slightly different from \(\kappa_1\) and \(\kappa_2\). Sketch the band structure for this case in the same figure. Is your sketch consistent with your answers to the previous two subquestions?

7. Suppose we set \(\kappa_3 = 0\). What set of simpler systems do we obtain? (it always helps to make a drawing!). What are the eigenfrequencies and their degeneracies of this set of simpler systems when \(\kappa_1=\kappa_2\)? (you can refer to the first exercise of the Bonds and Spectra lecture)