```python tags=["initialize"]
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)} .
$$

```python
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](https://en.wikipedia.org/wiki/Peierls_transition).

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:

```python
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](https://ocw.mit.edu/courses/physics/8-04-quantum-physics-i-spring-2013/study-materials/MIT8_04S13_OnCommEigenbas.pdf#page=3).

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