8_many_atoms_solutions.md

import matplotlib
from matplotlib import pyplot

import numpy as np

from common import draw_classic_axes, configure_plotting

configure_plotting()

pi = np.pi

Solutions for lecture 8 exercises

Quick warm-up exercises

The part first part of

$\omega^2$

will always be equal or larger than the second part. Therefore

$\omega^2 \geq 0$

.

Be cautious with the difference in the unit cell size.

Values of

$k$

outside of the 1st Brillouin zone describe the same solutions.

Exercise 1: analyzing the diatomic vibrating chain

Use the small angle approximation

$\sin(x) \approx x$

to ease calculations. For the Taylor polynomial take

$\omega^2 = f(x) \approx f(0) + f'(0)k + \frac{1}{2} f''(0)k^2$

(some terms vanish, computation is indeed quite tedious, but it's a 'fun' task). Calculating the group velocity yields

$$\left| v_g \right| = \sqrt{\frac{\kappa a^2}{2(m_1+m_2)}}$$

Optical branch corresponds with (+) in the equation given in the lecture notes. Be smart: you do not have to calculate the derivatives again. Finding the Taylor polynomial andcomputing the group velocity results in

$$\left| v_g \right|=0$$

Density of states is given as

$g(\omega) = dN/d\omega = dN/dk \times dk/d\omega$

. We know

$dN/dk = 2L/2\pi = L/\pi$

since we have 1D and positive and negative

$k$

-values.

$dk/d\omega$

can be computed using the group velocity:

$dk/d\omega = (d\omega/dk)^{-1} = (v_g)^{-1}$

Exercise 2: the Peierls transition

A unit cell consits of exactly one

$t_1$

and one

$t_2$

hopping.

Using the hint we find:

$$E \phi_n = \epsilon \phi_n + t_1 \psi_n + t_2 \psi_{n-1}$$

$$E \psi_n = t_1 \phi_n + t_2 \phi_{n+1} + \epsilon \psi_n$$

Notice that the hopping, in this case, is without the '-'-sign!

Using the Ansatz and rearranging the equation yields:

$$E \begin{pmatrix} \phi_0 \\ \psi_0 \end{pmatrix} = \begin{pmatrix} \epsilon & t_1 + t_2 e^{-ika} \\ t_1 + t_2 e^{ika} & \epsilon \end{pmatrix} \begin{pmatrix} \phi_0 \\ \psi_0 \end{pmatrix}$$

The dispersion is given by:

$$E = \epsilon \pm \sqrt{t_1^2 + t_2^2 + 2t_1t_2\cos(ka)} .$$

pyplot.figure()
k = np.linspace(-2*pi, 2*pi, 400)
t1 = 1;
t2 = 1.5;
pyplot.plot(k, -(t1+t2)*np.cos(k/2),'r',label='1 atom dispersion')
pyplot.plot(k[199:100:-1],-(t1+t2)*np.cos(k[0:99]/2),'r--',label='1 atom dispersion with folded Brillouin zone')
pyplot.plot(k[299:200:-1],-(t1+t2)*np.cos(k[300:399]/2),'r--')
pyplot.plot(k, np.sqrt(t1**2 + t2**2+2*t1*t2*np.cos(k)),'b',label='2 atom dispersion')
pyplot.plot(k, -np.sqrt(t1**2 + t2**2+2*t1*t2*np.cos(k)),'b')

pyplot.xlabel('$ka$'); pyplot.ylabel(r'$E-\epsilon$')
pyplot.xlim([-2*pi,2*pi])
pyplot.ylim([-1.1*(t1+t2),1.1*(t1+t2)])
pyplot.xticks([-2*pi, -pi, 0, pi,2*pi], [r'$-2\pi$',r'$-\pi$', 0, r'$\pi$',r'$2\pi$'])
pyplot.yticks([-t1-t2, -np.abs(t1-t2), 0, np.abs(t1-t2), t1+t2], [r'$-t_1-t_2$',r'$-|t_1-t_2|$', '0', r'$|t_1-t_2|$', r'$t_1+t_2$']);
pyplot.vlines([-pi, pi], -2*(t1+t2)*1.1,2*(t1+t2)*1.1, linestyles='dashed');
pyplot.hlines([-np.abs(t1-t2), np.abs(t1-t2)], -2*pi, 2*pi, linestyles='dashed');
pyplot.fill_between([-3*pi,3*pi], -np.abs(t1-t2), np.abs(t1-t2), color='red',alpha=0.2);

pyplot.legend(loc='lower center');

(Press the magic wand tool to enable the python code that created the figure to see what happends if you change

$t_1$

and

$t_2$

.)

Notice that the red shaded area is not a part of the Band structure anymore! Also check the wikipedia article.

For the group velocity we find

$$v_g(k) = \mp \frac{t_1 t_2}{\hbar}\frac{a \sin(ka)}{\sqrt{t_1^2 + t_2^2 + 2 t_1 t_2\cos(ka)}},$$

and for the effective mass we obtain

$$m^*(k) = \mp \frac{\hbar^2}{a^2 t_1 t_2}\left[ \frac{t_1 t_2 \sin^2(ka)}{(t_1^2 + t_2^2 +2 t_1 t_2\cos(ka))^{3/2}} + \frac{\cos(ka)}{\sqrt{t_1^2 + t_2^2 +2 t_1 t_2\cos(ka)}} \right]^{-1}.$$

We know

$g(E) = \frac{dN}{dk} \frac{dk}{dE} = \frac{L}{\pi} \frac{1}{\hbar v_g(E)}$

with

$v_g(k)$

from the previous subquestion. Rewriting

$v_g(k)$

to

$v_g(E)$

and substituting it in

$g(E)$

gives

$$g(E) = \frac{4L}{a \pi} \frac{\left| E-\epsilon \right|}{\sqrt{4t_1^2 t_2^2 -\left[ (E-\epsilon)^2 -t_1^2 - t_2^2 \right]^2}}.$$

Graphically the density of states looks accordingly:

pyplot.subplot(1,3,1)
k = np.linspace(-2*pi, 2*pi, 400)
t1 = 1;
t2 = 1.5;
pyplot.plot(k, -(t1+t2)*np.cos(k/2),'r',label='1 atom dispersion')
pyplot.plot(k[199:100:-1],-(t1+t2)*np.cos(k[0:99]/2),'r--',label='1 atom dispersion with folded Brillouin zone')
pyplot.plot(k[299:200:-1],-(t1+t2)*np.cos(k[300:399]/2),'r--')
pyplot.plot(k, np.sqrt(t1**2 + t2**2+2*t1*t2*np.cos(k)),'b',label='2 atom dispersion')
pyplot.plot(k, -np.sqrt(t1**2 + t2**2+2*t1*t2*np.cos(k)),'b')

pyplot.xlabel('$ka$'); pyplot.ylabel(r'$E-\epsilon$')
pyplot.xlim([-2*pi,2*pi])
pyplot.xticks([-2*pi, -pi, 0, pi,2*pi], [r'$-2\pi$',r'$-\pi$', 0, r'$\pi$',r'$2\pi$'])
pyplot.yticks([-t1-t2, -np.abs(t1-t2), 0, np.abs(t1-t2), t1+t2], [r'$-t_1-t_2$',r'$-|t_1-t_2|$', '0', r'$|t_1-t_2|$', r'$t_1+t_2$']);

pyplot.subplot(1,3,2)
w = np.sqrt(t1**2 + t2**2+2*t1*t2*np.cos(k))
pyplot.hist(w,30, orientation='horizontal',ec='black',color='b');
pyplot.hist(-w,30, orientation='horizontal',ec='black',color='b');
pyplot.xlabel(r'$g(E)$')
pyplot.ylabel(r'$E-\epsilon$')
pyplot.yticks([],[])
pyplot.xticks([],[])

pyplot.subplot(1,3,3)
w = -(t1+t2)*np.cos(k/2)
pyplot.hist(w,60, orientation='horizontal',ec='black',color='r');
pyplot.xlabel(r'$g(E)$')
pyplot.ylabel(r'$E-\epsilon$')
pyplot.yticks([],[])
pyplot.xticks([],[])

pyplot.suptitle('Density of states for 2 atom unit cell and 1 atom unit cell');

Exercise 3: atomic chain with 3 different spring constants

The unit cell should contain exactly one spring of

$\kappa_1$

,

$\kappa_2$

and

$\kappa_3$

and exactly three atoms.

The equations of motion are

\begin{align} m\ddot{u}{n,1} &= -\kappa_1(u{n,1} - u_{n,2}) - \kappa_3(u_{n,1} - u_{n-1,3})\ m\ddot{u}{n,2} &= -\kappa_2(u{n,2} - u_{n,3}) - \kappa_1(u_{n,2}-u_{n,1})\ m\ddot{u}{n,3} &= -\kappa_3(u{n,3} - u_{n+1,1}) - \kappa_2(u_{n,3}-u_{n,2}) \end{align}

Substitute the following ansatz in the equations of motion:

$$\begin{pmatrix} u_{1,n} \\ u_{2,n} \\ u_{3,n} \end{pmatrix} = e^{i\omega t - ikx_n} \begin{pmatrix} A_1 \\ A_2 \\ A_3 \end{pmatrix}$$

Because the spring matrix and the matrix

$X$

commute, they share a common set of eigenvectors. The eigenvalues of the matrix

$X$

are

$\lambda = -1$

with an eigenvector

$$\begin{pmatrix} 1 \\ 0 \\ -1 \end{pmatrix}$$

and

$\lambda = +1$

with eigenvectors

$$\begin{pmatrix} 1 \\ 0 \\ 1 \end{pmatrix}, \quad \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}.$$

These eigenvectors can be used to calculate the eigenvalues of the spring matrix. However, be cautious! The eigenvalue

$\lambda = +1$

is degenerate and to find the eigenvalue of the spring matrix we should take a linear combination of the two corresponding eigenvectors. This is explained at page 3 in this document.

The eigenvalues of the spring matrix are

$$\omega^2 = \begin{pmatrix} \omega_1^2 \\ \omega_2^2 \\ \omega_3^2 \end{pmatrix} = \frac{1}{m} \begin{pmatrix} q \\ k_3 + \frac{3q}{2} - \frac{\sqrt{4k_3^2 - 4k_3q + 9q^2}}{2} \\ k_3 + \frac{3q}{2} + \frac{\sqrt{4k_3^2 - 4k_3q + 9q^2}}{2} \end{pmatrix}$$

If

$\kappa_1 = \kappa_2 = \kappa_3$

then we have the uniform mono-atomic chain. If the length of the 3 spring constant unit cell is

$a$

, then the length of the mono-atomic chain is

$a/3$

.