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

Warm-up exercises

1. We rewrite the equation for \(\omega^2\) to the form of excercise 1. The term under the square root is always smaller than 1, from which we conclude that \(\omega^2 \geq 0\). This is important, because a negative value would yield imaginary solutions for \(\omega\), which does not correspond to an oscillating motion of the atoms (plug it into the Ansatz to see this)

2. Values of \(k\) that differ by an integer multiple of \(2\pi/a\) describe the same wave. Therefore, all information is already contained in the 1st Brillouin zone.

3. The LCAO wavefunction is \(|\Psi\rangle = \sum_n\phi_n|n\rangle\). The Schrödinger equation yields \(\varepsilon \phi_n = \varepsilon_0 \phi_n - t \phi_{n+1} - t \phi_{n-1}\). Plugging in the Ansatz \(\phi_n = \phi_0 e^{ikna}\) yields $$ \varepsilon = \varepsilon_0 -2t\cos(ka) $$

4. A band gap will open at \(k=\pi/a\), with \(a\) the size of the new unit cell (twice the interatomic distance)

Exercise 1: analyzing the diatomic vibrating chain

We have $$ \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), $$

1. When \(m_1=m_2=m\), we have \(\mu=m/2\). We get

   \[ \omega = \sqrt{\frac{2\kappa}{m}}\sqrt{\left(1\pm|\cos(ka/2)|\right)} =2\sqrt{\frac{\kappa}{m}} \begin{cases} |\sin(ka/4)| \\|\cos(ka/4)| \end{cases} \]

   We observe that this is identical to the dispersion of the monatomic chain derived in the tight-binding lecture, but folded back into the 1st Brillouin zone because the unit cell is twice as large.

2. The dispersion sketch should resemble that of the monatomic chain with identical masses, except for the band gap that opens at the edge of the Brillouin zone at \(ka=\pm\pi\) because of the different masses. The gap is centered at \(\sqrt{\tfrac{\kappa}{\mu}}\).

3. Because the \(\sin^2\) term in the square root is small, we can use the Taylor approximation \(\sqrt{1-x}\approx 1-x/2\). For the acoustic branch we get $$ \omega^2 \approx \frac{1}{2}\frac{\kappa}{\mu}\frac{4\mu^2}{m_1m_2}(ka/2)^2 $$ where we also used \(\sin(x) \approx x\). Therefore, the group velocity is $$ v_g=d\omega/ dk = \sqrt{\frac{\kappa a^2}{2(m_1+m_2)}} $$

4. For the optical branch, we get $$ \omega^2 \approx \frac{\kappa}{\mu}\left(2-\frac{1}{2}\frac{4\mu^2}{m_1m_2}(ka/2)^2\right) $$ from which we find that \(v_g(k=0) = d\omega/dk|_{k=0}=0\)

5. The density of states is \(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, and \(dk/d\omega\) is given by the inverse of the group velocity: \(dk/d\omega = (d\omega/dk)^{-1} = (v_g)^{-1}\) calculated in subquestion 2. We observe that, at small \(k\), the density of states is constant as in the 1D Debye model.

6. The band gap opens at \(ka=\pi\). When the masses are only slightly different, we can use \(\sqrt{1\pm x}\approx1\pm x/2\) to get $$ \omega_\pm=\sqrt{\frac{\kappa}{\mu}}\sqrt{1 \pm \sqrt{1-\frac{4\mu^2}{m_1m_2}}} \approx \sqrt{\frac{\kappa}{\mu}}\left(1 \pm \frac{1}{2}\sqrt{1-\frac{4\mu^2}{m_1m_2}}\right) $$ The band gap is \(\Delta\omega = \omega_+-\omega_-\)

Exercise 2: Analyzing the LCAO chain with alternating hoppings

1. The unit cell contains exactly two atoms. By looking at the figure, we can immediately write down the Schrodinger equation relating the probability amplitudes of the atomic orbitals: $$ E \phi_n = \varepsilon_0 \phi_n + t_1 \psi_n + t_2 \psi_{n-1} $$ $$ E \psi_n = \varepsilon_0 \psi_n + t_1 \phi_n + t_2 \phi_{n+1} $$

2. This follows directly from substituting the Ansatz

3. Solving the eigenvalue problem, we find the dispersion:

   \[ E_\pm = \epsilon \pm \sqrt{t_1^2 + t_2^2 + 2t_1t_2\cos(ka)} . \]

   This dispersion is similar to that of a chain with \(t_1=t_2\), except that band gaps form at the edge of the Brillouin zone (\(k=\pi/a\)). The band gap has size \(\Delta E = E_+-E_- = 2|t_1-t_2|\). The dispersion is plotted below

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

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

4. Realizing that the density of states is proportional to 1 over the slope of the dispersion (= the group velocity) enables a quick, graphical construction of the density of states.

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

To calculate the density of states, we use \(g(E) = \frac{dN}{dk} \frac{dk}{dE} = 2_s\frac{L}{\pi} \frac{1}{\hbar v_g(E)}\) with \(v_g(k)\). Rewriting \(v_g(k)\) to \(v_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');

6. The eigenvectors at \(k=0\) are \([\phi_0 \quad \psi_0] =\tfrac{1}{\sqrt{2}} [1 \quad \pm 1]\). The corresponding wavefunctions are \(|\Psi_\pm(k=0)\rangle = \tfrac{1}{\sqrt{2}}\sum_n(\phi_0|n,1\rangle \pm \psi_0|n,2\rangle)\). The (+) corresponds to the low-energy wavefuction. The wavefunctions are plotted below.

7. The eigenvectors at \(k=\pi/a\) are also \([\phi_0 \quad \psi_0] =\tfrac{1}{\sqrt{2}} [1 \quad \pm 1]\). The corresponding wavefunctions are \(|\Psi_\pm(k=\pi/a)\rangle = \tfrac{1}{\sqrt{2}}\sum_n(-1)^n(\phi_0|n,1\rangle \pm \psi_0|n,2\rangle)\). The wavefunctions are plotted below (dashed lines). We can see they are equal except for a shift. Therefore, if the hoppings are equal, they should be degenerate in energy. However, the unequal hoppings lift this degeneracy. The lowest energy is associated with the wavefunction that is an even function with respect to a bond with the largest hopping.

def gausspeak(x, mu=0, sigma=1, amp = 1):
    y = amp*np.exp(-(x-mu)**2/2/sigma**2)
    return y
N=400
x = np.linspace(-0.25,2.75,N)
groundstate = np.zeros(N)
excitedstate = np.zeros(N)
bandgap_lowerstate = np.zeros(N)
bandgap_upperstate = np.zeros(N)
muvals = np.linspace(0,2.5,6)
n=1
sigma = 0.07
for mu in muvals: 
  groundstate += gausspeak(x,mu=mu, sigma=sigma)
  excitedstate += (-1)**n*gausspeak(x,mu=mu, sigma=sigma)
  bandgap_lowerstate += (-1)**np.floor(n/2)*gausspeak(x,mu=mu, sigma=sigma)
  bandgap_upperstate += (-1)**np.floor((n+1)/2)*gausspeak(x,mu=mu, sigma=sigma)
  n+=1
pyplot.figure()
pyplot.plot(np.linspace(0,2.5,6),np.zeros(6),'ko')
pyplot.plot(x, groundstate,'b',label=r'$|\Psi_+(k=0)\rangle$')

pyplot.plot(np.linspace(0,2.5,6),np.zeros(6)+6,'ko')
pyplot.plot(x, excitedstate+6,'r',label=r'$|\Psi_-(k=0)\rangle$')

pyplot.plot(np.linspace(0,2.5,6),np.zeros(6)+2,'ko')
pyplot.plot(x, bandgap_lowerstate+2,'b--',label=r'$|\Psi_+(k=\pi/a)\rangle$')

pyplot.plot(np.linspace(0,2.5,6),np.zeros(6)+4,'ko')
pyplot.plot(x, bandgap_upperstate+4,'r--',label=r'$|\Psi_-(k=\pi/a)\rangle$')
pyplot.xlabel('$x/a$')
pyplot.ylim([-0.5, 7.5]);
pyplot.legend()
pyplot.suptitle('Wavefunctions at $k=0$ and $k=\pi/a$');
pyplot.yticks([],[]);

8. If we have one electron per atom, all states below \(E_F = \varepsilon_0\) are filled. Because a band gap appears when the hoppings become different, the total energy of the electrons will be lower than if the hoppings are equal. In contrast, for 2 electrons per atom, all electronic states would be filled with electrons and therefore the total energy would be insensitive to a difference in the hoppings.

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 into 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. Making use of the symmetry, we guess the form of the eigenvectors. As always, it is very useful to make a drawing. From a drawing of the unit cell, we expect an eigenvector \(\mathbf{v_1} = [1 \quad 1 \quad 1]^\text{T}\) (all atoms moving in the same direction) with eigenfrequency equal to zero, an eigenvector \(\mathbf{v_2} = [1 \quad a \quad 1]^\text{T}\) (outer atoms moving oppositely to the middle atom, with the middle atom having twice the amplitude), and an eigenvector \(\mathbf{v_3} = [1 \quad 0 \quad -1]^\text{T}\)

   We find the eigenvalues and the value of \(a\) by multiplying the spring matrix by these eigenvectors, yielding \(a=-2\) and

   \[ \omega^2 = \begin{pmatrix} \omega_1^2 \\ \omega_2^2 \\ \omega_3^2 \end{pmatrix} = \frac{\kappa}{m} \begin{pmatrix} 0 \\ 3 \\ 1+2\kappa_3/\kappa \end{pmatrix} \]

   where we defined \(\kappa = \kappa_1=\kappa_2\). More formally, we can use that 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 $ [1 \quad 0 \quad -1]^\text{T}$, and \(\lambda = +1\) with eigenvectors \([ 1 \quad 0 \quad 1]^\text{T}\) and \([0 \quad 1 \quad 0 ]^\text{T}\).

   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.

5. Trying the same \(\mathbf{v}_2\) and \(\mathbf{v_3}\) as in the previous subquestion, we find the eigenvalues

   \[ \omega^2 = \begin{pmatrix} \omega_1^2 \\ \omega_2^2 \\ \omega_3^2 \end{pmatrix} = \frac{\kappa}{m} \begin{pmatrix} 1 \\ \eta + \frac{3}{2} -\tfrac{1}{2} \sqrt{(1-2\eta)^2+8} \\ \eta + \frac{3}{2} +\tfrac{1}{2} \sqrt{(1-2\eta)^2+8} \end{pmatrix} \]

   where we defined \(\eta = \kappa_3/\kappa\)

6. If \(\kappa_1 = \kappa_2 = \kappa_3\) then we have the uniform mono-atomic chain. This mono-atomic chain has unit cell of size \(a/3\), so the Brillouin zone will span \(-3\pi/a\) to \(3\pi/a\). If we then make one of the spring constants slightly different, bandgaps will open at at \(k=pi/a\) and \(k=2\pi/a\), reflecting the new periodicity of the lattice.

7. When \(\kappa_3=0\), we obtain a set of isolated 3-atom molecules. We have calculated the eigenfrequencies and eigenvectors in exercise 1 of the Bonds and Spectra lecture.