<jupyter_start><jupyter_text>Electrons in a cristal - Bloch's theorem CrystalAn ideal crystal is made of periodic replicas of a structural unit, the **unit cell**. The periodicity is decribed as a **Bravais lattice**, *i.e.*, as multiples of **primitive vectors**:$$ \vec{R}_i = l_i \vec{a}_1 + m_i \vec{a}_2 + n_i \vec{a}_3$$where $l_i, m_i, n_i$ are integers ($l_i, m_i, n_i \in \mathbb{Z}^3$) and $\vec{a}_1, \vec{a}_2, \vec{a}_3$ are the primitive vectors of the unit cell. Reciprocal latticeFor several applications, such as diffraction problems, it can be extremly useful to work in the reciprocal space made of $\vec{K}$ vectors. For all $\vec{R}_i$ vectors of the real space lattice,$$e^{i\vec{K}\cdot\vec{R}_i} = 1$$It means a plane wave $e^{i\vec{K}\cdot\vec{r}}$ with the wavevector $\vec{K}$ has the same periodicity as the crystal. Here we will see that the reciprocal lattice is of key interest in the study of electronic structure of crystals.It is a Bravais lattice with primitive vectors $\vec{b}_i, i\in\{1,2,3\}$ verifying$$\vec{b}_i\cdot\vec{a}_j = 2\pi\delta_{ij}$$One can show that $\vec{b}_i$ vectors satisfying this property write as:$$\vec{b}_1 = 2\pi\dfrac{\vec{a}_2\wedge\vec{a}_3}{V}$$$$\vec{b}_2 = 2\pi\dfrac{\vec{a}_3\wedge\vec{a}_1}{V}$$$$\vec{b}_3 = 2\pi\dfrac{\vec{a}_1\wedge\vec{a}_2}{V}$$where $V$ is the volume of the unit cell: $V = \vec{a}_1\cdot(\vec{a}_2\wedge\vec{a}_3)$ Brillouin zoneThe Brillouin zone is a primitive cell of the reciprocal lattice defined as all the points closer to the $\Gamma$-point -- the origin $\vec{0}$ of the reciprocal lattice -- than from any other $\vec{K}$ vector.Its volume is $\dfrac{(2\pi)^3}{V}$. Bloch's theoremLet's consider eigenvectors $\psi(\vec{r})$ of the one-electron hamiltonian $\hat{H}$:$$\hat{H} = -\dfrac{\hbar^2}{2m}\vec{\nabla}^2 + U(\vec{r})$$where $U(\vec{r})$ is a potential with the same periodicity as the Bravais lattice:$$\forall \vec{R} \text{ of the Bravais lattice, }U(\vec{r}+\vec{R}) = U(\vec{r})$$Then, Bloch's theorem states that $\psi(\vec{r})$ writes:$$\psi(\vec{r}) = e^{i\vec{k}\cdot\vec{r}} u_{\vec{k}}(\vec{r})$$with $u_{\vec{k}}(\vec{r})$ a function with the same periodicity:$$\forall \vec{R} \text{ of the Bravais lattice, }u_{\vec{k}}(\vec{r}+\vec{R}) = u_{\vec{k}}(\vec{r})$$  ProofFor a lattice vector $\vec{R}$, we introduce the translation operator $\hat{T}_{\vec{R}}$:$$\hat{T}_{\vec{R}}\psi(\vec{r}) = \psi(\vec{r}+\vec{R})$$As $U(\vec{r})$ is periodic, the hamiltonian is translational-invariant. Then $\hat{H}$ and $\hat{T}_{\vec{R}}$ commute.**Q1/** Show that $\hat{T}_{\vec{R}}\hat{T}_{\vec{R}'} = \hat{T}_{\vec{R}+\vec{R}'}$ $\hat{T}_\vec{R}$ is a representation of the group of translations of vector $ \vec{R}$. As $\hat{H}$ and $\hat{T}_\vec{R}$ commute, it is possible to choose *common eigenvectors*:\begin{eqnarray}\hat{H}\psi(\vec{r}) = & \epsilon & \psi(\vec{r})\\\hat{T}_\vec{R}\psi(\vec{r}) = & c(\vec{R}) & \psi(\vec{r})\end{eqnarray} **Q2/** Show that the eigenvalues $c$ of the translation operator satisfy:$$c(\vec{R}+\vec{R}') = c(\vec{R})\times c(\vec{R}')$$ **Q3/** Show that a function $u_\vec{k}(\vec{r}) = e^{-i\vec{k}\cdot\vec{r}}\psi(\vec{r})$ has the same periodicity as the Bravais lattice. This achieves the proof of Bloch's theorem!  Additional remarks 1. $\vec{k}$ can always be chosen in the Brillouin zone.If $\vec{K}$ is any vector of the reciprocal lattice and $u(\vec{r}$ is a periodic function, then $e^{i\vec{K}\cdot\vec{r}}u(\vec{r})$ is periodic. So we can shift $\vec{k}$ towards the Brillouin zone using the appropriate $\vec{K}$ and writing:$$e^{i(\vec{k}+\vec{K})\cdot\vec{r}} u(\vec{r}) = e^{i\vec{k}\cdot\vec{r}}\times e^{i\vec{K}\cdot\vec{r}}u(\vec{r})$$2. We consider a Bloch state $\psi_\vec{k}(\vec{r}) = e^{i\vec{k}\cdot\vec{r}}u_\vec{k}(\vec{r})$, eigenvector for the hamiltonian $\hat{H}$:$$\hat{H}\psi_\vec{k}(\vec{r}) = \epsilon(\vec{k})\psi_\vec{k}(\vec{r})$$We can write,\begin{eqnarray}\hat{H}e^{i\vec{k}\cdot\vec{r}}u_\vec{k}(\vec{r}) = & -\dfrac{\hbar^2}{2m}\vec{\nabla}^2e^{i\vec{k}\cdot\vec{r}}u_\vec{k}(\vec{r}) + U(\vec{r})e^{i\vec{k}\cdot\vec{r}}u_\vec{k}(\vec{r})\\= & e^{i\vec{k}\cdot\vec{r}} \left[ \dfrac{\hbar^2}{2m}(-i\vec{\nabla} + \vec{k})^2 + U(\vec{r}) \right]u_\vec{k}(\vec{r})\end{eqnarray}Then, $u_\vec{k}(\vec{r})$ is an eigenvector of $\hat{H}_\vec{k}$ with the eigenvalue $\epsilon(\vec{k})$ where$$\hat{H}_\vec{k} = \dfrac{\hbar^2}{2m}(-i\vec{\nabla} + \vec{k})^2 + U(\vec{r})$$Since $u_\vec{k}(\vec{r})$ is periodic, the hamiltonian $\hat{H}_\vec{k}$ has usually discrete eigenvalues. We use the integer $n$ to order them:$$\hat{H}_\vec{k}u_{n,\vec{k}}(\vec{r}) = \epsilon_n(\vec{k})u_{n,\vec{k}}(\vec{r})$$which is the same as$$\hat{H}\psi_{n,\vec{k}}(\vec{r}) = \epsilon_n(\vec{k})\psi_{n,\vec{k}}(\vec{r})$$with $$\psi_{n,\vec{k}} = e^{i\vec{k}\cdot\vec{r}}u_{n,\vec{k}}(\vec{r})$$3. Although $\vec{k}$ can always be chosen within the Brillouin zone, the eigenvalues $\epsilon_n(\vec{k})$ are usually periodically replicated on the whole reciprocal space:For any reciprocal lattice vector $\vec{K}$,$$\epsilon_n(\vec{k}+\vec{K}) = \epsilon_n(\vec{k})$$4. As we expect $\epsilon_n(\vec{k})$ to be continuous with $\vec{k}$, then $\epsilon_n(\vec{k})$ necessarily admits a minimum and a maximum. We can then talk of an **energy band**.  1D application to polyacetylene ModellingWe consider an ideal 1D polymere such as polyacetylene which is made of $N$ periodic replicas of the $(CH_2=CH_2-)$ unit. To make it a general 1D model, we will name A and B the groups of atoms in the unit cell. Thus the model can be polyacetylene but also boron nitride BN or any kind of substituted polyacetylene, $ (CH_2=CCl_2-)$ for example. We will only consider the $p_z$ electrons of A and B. This means they are $2$ active electron per unit cell, 1 per group (A and B). We use $d$ for the length of the unit cell in real space. Interactions between atom groupsWe name $\left\vert A,n\right>$ and $\left\vert B,n\right>$ the $p_z$ orbital on respectively A and B in the $n^\mathrm{th}$ cell. We use in this exemple the tight-binding model which is similar to the Hückel model. The groups A and B interact only with neighboring groups:$$\left = \alpha_1\delta_{n,m}$$$$\left = \alpha_2\delta_{n,m}$$$$\left = \beta_1\delta_{n,m}$$$$\left = \beta_2\delta_{n,m-1}$$where $\alpha_1$ and $\alpha_2$ are the energies of the $p_z$ orbital of A and B (resp.) and $\beta_1$ and $\beta_2$ quantify the coupling between neighboring A and B groups. Brillouin zone **Q4/** What is the length of the unit cell in reciprocal space? What is a vector $\vec{k}$ in the Brillouin zone?  Cristal orbitalsAn electronic orbital for this system reads as a linear combination of the $p_z$ orbitals:$$\left\vert\Psi\right> = \dfrac{1}{\sqrt{N}}\sum_n\left( c_n\left\vert A,n\right> + d_n\left\vert B,n\right>\right)$$ **Q5/** As $\left\vert\Psi\right>$ is normalized, write the condition on the $c_n$ and $d_n$ coefficients. After Bloch's theorem, the eigenstates of the hamiltonien $\hat{H}$ of this system reads$$\left\vert\Psi_k\right> = \dfrac{1}{\sqrt{N}}\sum_n e^{ikdn}c(k)\left\vert A,n\right> + e^{ikd(n+\frac{1}{2})}d(k)\left\vert B,n\right>$$where $-\dfrac{\pi}{d}<k<\dfrac{\pi}{d}$. The periodic orbital is:$$\left\vert u_k\right> = \dfrac{1}{\sqrt{N}}\sum_n c(k)\left\vert A,n\right> + d(k)\left\vert B,n\right>$$NB: This definition means that the group B is in the middle of the cell (at $d/2$ from A). Any other translation is equivalent, they are simply different definitions of $c(k)$ and $d(k)$. Solution of the time-independent Schrödinger equation$\left\vert\Psi_k\right>$ is solution of the time-independent Schrödinger equation:$$\hat{H}\left\vert\Psi_k\right> = \epsilon(k)\left\vert\Psi_k\right>$$ **Q6/** Using the model described previously, compute $\hat{H}\left\vert A,n\right>$ and $\hat{H}\left\vert B,n\right>$. **Q7/** Simplify $\hat{H}\left\vert\Psi_k\right>$. It follows that:$$\hat{H}_k\left\vert u_k\right> = \dfrac{1}{\sqrt{N}}\sum_n \left( \alpha_1 c(k) + e^{ikd/2} \beta_1 d(k) + e^{-ikd/2} \beta_2 d(k)\right) \left\vert A,n \right> + \left( e^{-ikd/2} \beta_1 c(k) + e^{ikd/2} \beta_2 c(k) + \alpha_2 d(k) \right) \left\vert B,n \right> = \epsilon(k)\left\vert u_k\right>$$We choose to represent $\left\vert u_k\right>$ as the vector $\begin{pmatrix} c(k)\\d(k)\end{pmatrix}$ which verifies the equation in "$n$-representation"$$H_k\begin{pmatrix} c(k)\\d(k)\end{pmatrix} = \begin{pmatrix} \alpha_1 & e^{ikd/2}\beta_1 + e^{-ikd/2}\beta_2\\ e^{-ikd/2}\beta_1 + e^{ikd/2}\beta_2 & \alpha_2\end{pmatrix}\begin{pmatrix} c(k)\\d(k)\end{pmatrix} = \epsilon(k)\begin{pmatrix} c(k)\\d(k)\end{pmatrix}\tag{E1}$$ Eigenvalues and band structureAn exact solution for the eigenvalues of the $2\times 2$ matrix can be found:\begin{eqnarray}&\epsilon_+(k) & = & \frac{\alpha_1 + \alpha_2}{2} + \sqrt{\frac{(\alpha_1-\alpha_2)^2}{4}+(\beta_1^2+\beta_2^2+2 \beta_1\beta_2 \cos kd)}\\\mathrm{and}&\\&\epsilon_-(k) & = & \frac{\alpha_1 + \alpha_2}{2} - \sqrt{\frac{(\alpha_1-\alpha_2)^2}{4}+(\beta_1^2+\beta_2^2+2 \beta_1\beta_2 \cos kd)}\end{eqnarray}This means that in general the electrons of the system are in two bands $\epsilon_+(k)$ and $\epsilon_-(k)$ which are continuous functions of $k$.  Numerical application<jupyter_code># required modules

import numpy as np
import matplotlib.pyplot as plt
from wannier_toolbox import *<jupyter_output><empty_output><jupyter_text>**Q8/** Define the $H_k$ 2x2 matrix by filling the following function:<jupyter_code>def hamiltonian(alpha1,alpha2,beta1,beta2,d,k):
    ham = np.zeros((2,2),dtype=complex)
    ...
    return ham
# parameters for numerical calculations

d = 2
alpha1, alpha2 = -1.0, -1.0
beta1, beta2 = -1.0, -1.0

N = 20 # number of k-points
# k-space
kpt = np.linspace(-np.pi/d, np.pi/d, N+1)
# exact solution

def eplus(k,alpha1,alpha2,beta1,beta2,d):
    return (alpha1+alpha2)/2 + np.sqrt((alpha1-alpha2)**2/4+(beta1**2+beta2**2+2*beta1*beta2*np.cos(k*d)))

def eminus(k,alpha1,alpha2,beta1,beta2,d):
    return (alpha1+alpha2)/2 - np.sqrt((alpha1-alpha2)**2/4+(beta1**2+beta2**2+2*beta1*beta2*np.cos(k*d)))<jupyter_output><empty_output><jupyter_text>**Q9/** Solve numerically the equation (E1). **Q10/** Plot and compare the exact and numerical solutions.  Periodic boundary conditionsWe now have a good idea the energy levels for this system. The next logical step is to find out whether they are occupied or not.At $T=0~\mathrm{K}$ and with the help of Fermi-Dirac statistics, one can say that levels with energies smaller than Fermi energy are occupied by $2s + 1$ electrons where $s$ is the spin ($s = \frac{1}{2}$). Thus, the Fermi energy only depends on the total number of electrons in the cristal. Similarly to the ideal electron gas, we therefore need to consider the boundary conditions.Let's consider a cristal with $N_1$ unit cells along $\vec{a}_1$, $N_2$ along $\vec{a}_2$ and $N_3$ along $\vec{a}_3$. This system has $N=N_1\times N_2\times N_3$ cells and a total volume $V_{tot} = N\times V$ where $V$ is the volume of the unit cell.Periodic boundary conditions state that:$$\forall i\in\{ 1,2,3\}, \Psi(\vec{r}+N_i\vec{a}_i) = \Psi(\vec{r})$$ **Q11/** Show that it implies that $\vec{k}$ reads$$\vec{k} = \sum_{i=1}^3\dfrac{m_i}{N_i}\vec{b}_i$$with $(m_1, m_2, m_3) \in \mathbb{Z}^3$. The reciprocal space unit cell volume can then be computed:\begin{eqnarray}\Delta^3\vec{k} & = & \dfrac{\vec{b}_1}{N_1}\cdot(\dfrac{\vec{b}_2}{N_2}\wedge\dfrac{\vec{b}_3}{N_3}\\& = & \dfrac{1}{N}\vec{b}_1\cdot (\vec{b}_2\wedge\vec{b}_3)\\& = & \dfrac{(2\pi)^3}{NV} = \dfrac{(2\pi)^3}{V_{tot}}\end{eqnarray}In the thermodynamic limit, we can replace the sum on $\vec{k}$ for an integral:$$\sum_{\vec{k}}\rightarrow\int\mathrm{d}^3\vec{k}\dfrac{V_{tot}}{(2\pi)^3}$$With $f(\epsilon)$ the Fermi function, the number of electrons in band $n$ is$$N_n = (2s + 1)\int\mathrm{d}^3\vec{k}\dfrac{V_{tot}}{(2\pi)^3}f(\epsilon_n(\vec{k}))$$And the number of electrons per cell is$$\dfrac{N_n}{N} = (2s+1)\dfrac{V}{(2\pi)^3}\int\mathrm{d}^3\vec{k}f(\epsilon_n(\vec{k}))$$The maximal number of electrons per cell can be find if $\forall n,\vec{k}, f(\epsilon_n(\vec{k})) = 1$ which means that all bands are fully occupied. This gives$$\dfrac{N_n}{N} = 2$$Similarly to molecular orbitals, bands are filled with a maximum of $2$ electrons per cell.As for polyacetylene, the band $\epsilon_-(\vec{k})$ is fully occupied and the band $\epsilon_+(\vec{k})$ is empty. **Q12/** Play with the model parameters $\alpha_1,\alpha_2,\beta_1,\beta_2$ to understand and illustrate the ideas of insulation and conduction in metals  Bloch states: phase and representationThe $2\times 2$ eigenvectors matrix of $H_k$ gives $c(k)$ and $d(k)$ for each band (one band per column) which is a representation of $\left\vert u_{n,k}\right>$. The phase of these complex numbers can be freely chosen as long as it verifies the condition $\vert c(k)\vert^2+\vert d(k)\vert^2 = 1$:$$c(k) \rightarrow e^{i\varphi(k)}c(k) \qquad \text{and} \qquad d(k) \rightarrow e^{i\varphi(k)}d(k)$$Since $H_k\in\mathcal{M}_{2,2}^\mathbb{C}$, it may not be possible to find a $\varphi(k)$ such as both $c(k)$ and $d(k)$ are real numbers. In the following we choose $\varphi(k)$ so $c(k)\in\mathbb{R}$.<jupyter_code>N=10

ck=np.zeros((N+1,2),dtype=complex)
dk=np.zeros((N+1,2),dtype=complex)
eps=np.zeros((N+1,2))
phi=np.zeros((N+1,2))
kpt=np.linspace(-np.pi/2,np.pi/2,N+1)
for ik, k in enumerate(kpt):
    H = hamiltonian(alpha1,alpha2,beta1,beta2,d,k)
    e,v = np.linalg.eigh(H)
    eps[ik,:] = e
    phi[ik,:] = np.angle(v[0,:]) # angle of the complex number c(k) in the complex plane
    ck[ik,:] = v[0,:]*np.exp(-1j*phi[ik,:])
    dk[ik,:] = v[1,:]*np.exp(-1j*phi[ik,:])<jupyter_output><empty_output><jupyter_text>The true Bloch states are not really the $\left\vert u_{n,k}\right>$ but rather the wavefunctions $\left\vert\Psi_{n,k}\right>$ as shown previously which has a periodicity of $N$ cells so it will be represented on these cells:<jupyter_code>x = np.linspace(0,(2*N-1)*d/2,2*N) # positions of A and B

ik = 6 # select the wavevector
print(f"For a k-point k = {kpt[ik]/np.pi:.2f}π")

band = 0 # choose the band, 0 = lower energy, 1 = higher energy

Psi = np.zeros(2*N,dtype=complex)

for cell in range(N):
    Psi[2*cell] = ck[ik,band]*np.exp(1j*kpt[ik]*x[2*cell])
    Psi[2*cell+1] = dk[ik,band]*np.exp(1j*kpt[ik]*x[2*cell+1])

plt.bar(x,np.real(Psi))
plt.xlabel("x")
plt.ylabel("Re(Psi)")
plt.show()

plt.bar(x,np.imag(Psi))
plt.xlabel("x")
plt.ylabel("Im(Psi)")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20353013614284737, 'fraction_numerical': 0.02079770130669768, 'mean_word_length': 4.9062626262626265, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '49426', 'license': 'CC0-1.0', 'n_tokens_mistral': 5635, 'n_tokens_neox': 5362, 'n_words': 1504, 'path': 'Bloch_Wannier.ipynb', 'repo_name': 'axgmz/Bloch_theorem_Notebook_2ADENS', 'size': '14617'}