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

1.

The part first part of \(\omega^2\) will always be equal or larger than the second part. Therefore \(\omega^2 \geq 0\).

2.

Be cautious with the difference in the unit cell size.

3.

Values of \(k\) outside of the 1st Brillouin zone describe the same solutions.

Exercise 1: analyzing the diatomic vibrating chain

1.

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)}} \]

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 \]

3.

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

1.

A unit cell consits of exactly one \(t_1\) and one \(t_2\) hopping.

2.

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!

3.

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} \]

4.

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.

5.

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}. \]

6.

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

1.

The unit cell should contain exactly one spring of \(\kappa_1\), \(\kappa_2\) and \(\kappa_3\) and exactly three atoms.

2.

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}\]

3.

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} \]

4.

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} \]

5.

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\).