    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)

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.

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} = - κ (u_{n} - u_{n - 1}) - κ (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 = N a$ , we have $u_{N} = u_{0}$ .

Electrons¶

The following figure shows the interatomic potential with atoms placed at $x_{k} = k a$ with $k \in 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)

png

Similarly to the triatomic system case, we formulate the molecular orbital via the LCAO model: $| Ψ ⟩ = \sum_{n} ϕ_{n} | n ⟩ .$ We assume only nearest-neighbor hopping $- t$ and an on-site energy $E_{0}$ . The coupled Schrödinger equation of the $| n ⟩$ orbital is: $E ϕ_{n} = E_{0} ϕ_{n} - t ϕ_{n + 1} - t ϕ_{n - 1} .$

Again, the periodic boundary conditions imply $ϕ_{N} = ϕ_{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} = A e^{i ω t - i k x_{n}},

and the ansatz for electrons is given similarly by

ϕ_{n} = B e^{i E t / ℏ - i k x_{n}} .

We wrote $x_{n} = n a$ 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 π / L$ . In the case of the electron system, periodic boundary conditions give $ϕ_{0} = ϕ_{N}$ , which results in

1 = e^{i k 0} = e^{i k N a} .

The above equation defines the allowed values of $k$ :

k = \frac{2 π p}{N a}, with p \in Z .

We use the quantized values of $k$ in our plane wave ansatz: $e^{i k x_{n}} = e^{i p \frac{2 π}{N a} n a} = e^{i \frac{2 π n p}{N}}$ . We notice something interesting if we investigate the case of $p \to p + N$ . In this case, the plane wave ansatz becomes $e^{i \frac{2 π n (p + N)}{N}} = e^{i \frac{2 π n p}{N} + i 2 π n} = e^{i \frac{2 π n p}{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');

png

How many different solutions did we expect to find?

We have a system with $N$ degrees of freedom (either $u_{n}$ or $ϕ_{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 ω^{2} A e^{i ω t - i k x_{n}} = κ A e^{i ω t} (- 2 e^{- i k x_{n}} + e^{- i k x_{n} + i k a} + e^{- i k x_{n} - i k a}) .

Searching for solutions with $A \neq 0$ we obtain

- m ω^{2} = κ (- 2 + e^{i k a} + e^{- i k a}) = κ [- 2 + 2 \cos (k a)] .

Or after a further simplification:

ω = \sqrt{\frac{2 κ}{m}} \sqrt{1 - \cos (k a)} = 2 \sqrt{\frac{κ}{m}} | \sin (k a / 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);

png

The periodicity of the dispersion relation is a consequence of what we observed before: since plane waves with $k$ -vectors differing by $2 π / a$ are exactly the same, the repeated periods of the dispersion relation describe the same plane waves.

Comparison to the Debye model¶

Sound velocity: At small $k$ , $\sin (k a / 2) \approx k a / 2$ . Therefore $ω \approx \sqrt{κ / m} k a = v_{s} k$ , with $v_{s}$ the sound velocity. We therefore justify the linear dispersion approximation in the Debye model!

Cut-off frequency: The Debye model introduces the cutoff frequency $ω_{D}$ to limit the number of phonon modes, and does not identify the origin of the cutoff. Now because of the finite number of plane waves, the integration over the $k$ -space has a finite size: $\sum_{p} \to \frac{L}{2 π} \int_{- π / a}^{π / a} d k$ . This automatically gives us a maximal frequency without additional assumptions.

Electrons¶

Once again we substitute the plane wave ansatz into the equations of motion:

E A e^{i E t / ℏ - i k n a} = E_{0} A e^{i E t / ℏ - i k n a} - t A e^{i E t / ℏ - i k (n + 1) a} - t A e^{i E t / ℏ - i k (n - 1) a},

Again, we are not interested in a trivial solution, hence we assume $A \neq 0$ and thus

E = E_{0} - t e^{- i k a} - t e^{i k a} = E_{0} - 2 t \cos (k a),

which gives us the dispersion relation below.

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$']);

png

We see that the electron dispersion consists of a band of allowed energies $E_{0} - 2 t < E < E_{0} + 2 t$ . That particles occupy bands of allowed energies is why a dispersion relation is also often called a band structure. Due to the spin degeneracy, each band has 2 $N$ possible states if we consider a system with $N$ atoms.

If each atom contains 2 electrons and a single orbital, all the states in the band must be occupied by electrons. Because all the available states are occupied, there is always exactly the same number of electrons moving in the positive direction, as there are in the negative. Hence, no matter what we do, our system is incapable of conducting electrons, and therefore we have derived the existence of insulators!

Let us also compare the electron band structure with the free electron model. Focusing on the dispersion relation close to the band bottom at $k = 0$ , we approximate the energy as $E \approx E_{0} - 2 t + t (k a)^{2} .$

If we compare this to the dispersion relation $E = ℏ k^{2} / 2 m$ of the free electrons model, we see that the band structure is similar, but the lowest available energy is $E_{0} - 2 t$ instead of $0$ , and the electrons behave as if they had a different mass $m^{*} = ℏ^{2} / 2 t a^{2}$ .

Group velocity, effective mass, density of states¶

(here we only discuss electrons; for phonons everything is the same except for replacing $E = ℏ ω$ )

Let us think what happens if we apply an external electric field to the crystal:

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) + .2 * (x - 1.5))
ax.plot(x, .2 * (x - 0.25), '--')
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)

png

The full Hamiltonian of the system is

H = \frac{p^{2}}{2 m} + U_{atomic} (x) + e E x,

where $U_{atomic}$ is the potential created by the nuclei, and $E$ the electric field.

A typical electric field is much smaller than the interatomic potential, and therefore we can start by obtaining the dispersion relation $E (k)$ without electric field (by applying the LCAO method), and then solve $H = E (k) + e E x .$

To derive how particles with an arbitrary dispersion relation move, we recall the Hamilton's equations for particle velocity $v$ and force $F$ :

\begin{aligned} v \equiv \frac{d r}{d t} & = \frac{\partial H (p, r)}{\partial p} \\ F \equiv \frac{d p}{d t} & = - \frac{\partial H (p, r)}{\partial r} . \end{aligned}

Substituting $p = ℏ k$ into the first equation we arrive to the expression for the electron group velocity $v \equiv ℏ^{- 1} \partial E / \partial k$ . From the second equation we obtain that the force acting on electron in a band stays $- e E$ , which in turn gives results in the acceleration

\frac{d v}{d t} = \frac{\partial v}{\partial p} \frac{d p}{d t} = F / m .

Comparing this expression with $d v / d t = F / m$ , we arrive to the effective mass:

m^{*} \equiv {(\frac{\partial v}{\partial p})}^{- 1} = {(\frac{\partial^{2} E}{\partial p^{2}})}^{- 1} = ℏ^{2} {(\frac{\partial^{2} E}{\partial k^{2}})}^{- 1} .

The group velocity describes how quickly electrons with a certain $k$ -vector move, while the effective mass describes how hard they are to accelerate by applying external force.

By using the dispersion relation we derived earlier, we obtain the effective mass like this:

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^*$')
pyplot.xticks([-pi, 0, pi], [r'$-\pi$', 0, r'$\pi$']);

png

Notice that the effective mass can be negative, which implies the electrons accelerate in the direction opposite to the applied force.

Density of states¶

The DOS is the number of states per unit energy. In 1D we have

g (E) = \frac{L}{2 π} \sum | d k / d E | = \frac{L}{2 π} \sum | v |^{- 1}

The sum goes over all possible values of $k$ and spin which have the same energy $E$ . If we are working in two or more dimensions, we must integrate over the values of $k$ with the same energy. Also take note that for energies below $E_{0} - 2 t$ or above $E_{0} + 2 t$ , there are no values of $k$ with that energy, so there is nothing to sum over.

Once again, starting from

E = E_{0} - 2 t \cos (k a),

we get

k a = \pm \arccos [(E - E_{0}) / 2 t],

and

| v |^{- 1} = | \frac{d k}{d E} | = \frac{1}{a} \frac{1}{\sqrt{4 t^{2} - (E - E_{0})^{2}}} .

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

We now add together the contributions of the positive and the negative momenta as well both spin orientations, and arrive to the density of states

g (E) = \frac{L}{2 π} \frac{4}{a} \frac{1}{\sqrt{4 t^{2} - (E - E_{0})^{2}}} .

A quick check: when the energy is close to the bottom of the band, $E = E_{0} - 2 t + δ E$ , we get $g (E) \propto δ E^{- 1 / 2}$ , as we expect in 1D.

The process of calculating the DOS at a given energy $E$ of a spin-independent Hamiltonian is done systematically with the following steps:

At a given energy $E$ , determine all of the values of $k$ which correspond to that $E$ using the dispersion relation.
Compute $| d k / d E |$ . Do this either by writing $k$ as a (multi-valued) function of $E$ and differentiating, or by computing $(d E / d k)^{- 1}$ .
Sum or integrate $d k / d E$ over the allowed values of $k$ found in 1 and multiply by any degeneracies (spin/polarization).
Multiply by spin degeneracy.

If the Hamiltonian depends on spin, then there is no spin degeneracy and the spin number $s$ must be treated in the same way as $k$ .

Summary¶

By using plane waves in real space as an Ansatz we found all normal modes and energies for phonons and electrons in 1D
Computing dispersion relations explains 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.

Exercises¶

Warm-up questions¶

Compare the expression of the effective mass with Newton's second law. Do you observe any similarities?
Check the units of the group velocity. Is it what you expect?
Do the same for the effective mass.
Calculate the effective mass of the free-electron dispersion relation. Is this what was expected?
Under what condition is the effective mass the same for each electron?

Exercise 1: Lattice vibrations¶

From the dispersion relation of a 1D monatomic chain given in the lecture notes, calculate the group velocity $v_{g}$ .
Using the group veolcity found above, calculate the density of states $g (ω)$ .
Sketch them.
From the 1D dispersion relation $ω (k)$ in the picture below, sketch the group velocity $v_{g} (k)$ and the density of states $g (ω)$ .

Exercise 2: Vibrational heat capacity of a 1D monatomic chain¶

Give the total energy $U$ of a 1D monatomic chain as an integral expression. To do so, first derive the density of states from the phonon dispersion relation derived in the lecture notes.
Give an integral expression for the heat capacity $C$ .
Compute the heat capacity numerically, using e.g. Python.
Do the same for $C$ in the Debye model and compare the two. What differences do you see?

Exercise 3: Next-nearest neighbors chain¶

Consider electrons in a 1D atomic chain again:

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)

png

Let's expand the one-dimensional chain model by extending the range of the interaction further than the nearest neighbors:

⟨ ϕ_{n} | H | ϕ_{n + 2} ⟩ \equiv - t^{'} \neq 0.

Write down the new Schrödinger equation for this system.

hint

There are now two more terms in the equation: $- t^{'} ϕ_{n - 2} - t^{'} ϕ_{n + 2}$ .
Solve the Schrödinger equation to find the dispersion relation $E (k)$ .

hint

Use the same Ansatz as for the nearest neighbors case: $ \phi_n = \phi_0 \exp(ikna) $.
Calculate the effective mass $m^{*}$ .
Sketch the effective mass as a function of $k$ for the cases $t = 2 t^{'}$ , $t = 4 t^{'}$ , and $t = 10 t^{'}$ .