In [1]:
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

# from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
%matplotlib widget
plt.rcParams['figure.figsize'] = [8,6]
plt.rcParams['figure.dpi'] = 100 # 200 e.g. is really fine, but slower

The emergence of band structures¶

If you have questions about this notebook, e-mail Oleg Malanyuk (oleg.malanyuk@epfl.ch)

This jupyter notebook is intended to provide a illustration of the band structures in solides, starting from the points of view of free electrons and atomic orbitals. The goals are:

  • To understand how the Bloch theorem leads to band structures.
  • To see how the electronic bands emerge from the interplay of periodic potential and kinetic energy.
  • Understanding the relationship between the Schrodinger equation and the calculation*

The Bloch theorem¶

In solids, the lattice has discrete translational symmetry. The Bloch theorem says that, such symmetry applies to the wavefunction as well, such that $$[H,R_i]=0$$ To study the energy levels in solides, one can adapt the simplest Hamiltonian:$$ H = -\frac{h^2}{2m} \nabla^2 + U(r)$$ where $U(r+R_i)=U(r)$

From the commutation relations of the Hamiltonian, the wavefunctions should be of the form $\psi(nR_i + r) = \psi(r) \exp(-inR_i)$ which is called Bloch functions.

From the atomic orbitals/free electrons to band structures¶

The Hamiltonian $$ H = -\frac{h^2}{2m} \nabla^2 + U(r)$$ contains two components: the kinetic energy $E_k=-\frac{h^2}{2m} \nabla^2$ and the potential $U$. Thus, we can look into the formation of band structure from the limits of free electrons($E_k >> U$) or localized orbitals($E_k<<U$) and treating the other term as a perturbation.

In [2]:
#setting up the starting parameters to initialise the plot
t_g = 0.5 #weight of kintetic energy
u_g = 0.5 #weight of potential energy
U_type_g = 'Sine Wave' #time of the potential
n_G=20 #number of G vectors, utilised in the size of the Hamiltonian matrix
#(technically should be infinite, but as we are concerned with only a few initial bands finite size is fine)
n_k=30 #number of k vectors, technically depends on number of unit cells N, the size of the system 

#the function below deals with changing potential type
def interactive_bands_V(U_type):
    '''
    :param k: wave number
    '''
    #Taking Fourier transform of the potential takes quite a bit of computational time, hence the interactive
    #part is split between changing the potential type and changing weights. To keep variables consistent
    #between the two functions, global variables are used
    global t_g, u_g, U_type_g, n_G, n_k 
    global ax1, ax2, ax3
    U_type_g = U_type #as potential type is changed, update the global variable
    t = t_g #read the weights of KE and PE from global variables
    u = u_g
    mat=np.zeros((n_G+1,n_G+1)) #predefine the size of the Hamiltonian matrix
    
    n = 30 #these two variables are used in plotting the potential and its Fourier transform
    fs = 500 #n is number of atoms (system size) and fs is number of points to plot between each atom
    x = np.linspace(-n, n, n*fs+1, endpoint=True) #position vector to be untilised in plotting the potential
    
    k_values=np.zeros((n_k+1,n_G+1)) #predefine the size of variables which will be used for band plots
    Energy_Eigenvalues=np.zeros((n_k+1,n_G+1))
    
   
    
    if (U_type == 'Sine Wave'):
        V = -np.cos(x*2*np.pi) #define potential, in this case cos()

        for i in range(n_k+1):
            mat=np.zeros((n_G+1,n_G+1))
            U=np.ones(n_G+1)*(-u) #define fourier components of the potential. As potential is cos() hence only G'-G = +-1 compontents will be present
            K=np.ones(n_G+1) #define dimensions of kinetic energy vector
            for j in range(n_G+1):
                K[j] = t*(i/n_k - 1/2 + j - n_G/2)**2 #for every
            mat=mat+np.diag(U,-1)[:-1,:-1]+np.diag(U,1)[:-1,:-1]+np.diag(K,0) 
            #put the kinetic terms on the diagonal of the Hamiltonian, and add the potential as off diagonal elements
            Energy_Eigenvalues[i]=np.linalg.eigvalsh(mat) #energy levels will be the eigenvalues of hamiltonian
            k_values[i]=np.ones(n_G+1)*(i-n_k/2)/n_k #keep track of which k value is the eigenvalues for

    elif (U_type == 'Series of Delta Functions'):
        V = (np.ones(x.shape) - np.ceil(x-np.floor(np.round(x,4)))) #convinient way to define series of delta functions
        
        for i in range(n_k+1):
            K=np.ones(n_G+1)
            for j in range(n_G+1):
                K[j] = t*(i/n_k - 1/2 + j - n_G/2)**2
            mat=np.ones((n_G+1,n_G+1))*(u) #in this case, all G'-G values are u, hence we just define all matrix elements to be u
            mat=mat+np.diag(K,0) #add kinetic term to the diagonal
            Energy_Eigenvalues[i]=np.linalg.eigvalsh(mat) #energy levels will be the eigenvalues of hamiltonian
            k_values[i]=np.ones(n_G+1)*(i-n_k/2)/n_k 
            
    elif (U_type == 'Square Wave'):
        V = 2*np.round((-np.cos(x*2*np.pi)+1)/2, 0) - np.ones(x.shape) #convinient way to define square wave
        
        for i in range(n_G+1):
            for j in range(n_G+1):
                if ((i - j) %2!=0):
                    mat[i][j] = (-1)**(int((i-j+1)/2))*u/(i-j) 
        #in this case, the fourier coefficients will be non-zero only if G'-G is odd, are proportional to 1/(G'-G)
        #and vary in sign, -1 if (G'-G)/2pi = = 1 (mod 4), 1 if (G'-G)/2pi = -1 (mod 4) [a=1 for convinience]
        for i in range(n_k+1):
            mat_sq = np.copy(mat)
            K=np.ones(n_G+1)
            for j in range(n_G+1):
                K[j] = t*(i/n_k - 1/2 + j - n_G/2)**2
            mat_sq=mat_sq+np.diag(K,0) #add kinetic energy to the hamiltonian
            Energy_Eigenvalues[i]=np.linalg.eigvalsh(mat_sq) #energy levels will be the eigenvalues of hamiltonian
            k_values[i]=np.ones(n_G+1)*(i-n_k/2)/n_k 
    
    #find fourier transform of the potential
    Vq = np.fft.fft(V)
    k = np.fft.fftfreq(n = x.shape[0], d = 2/fs)
    
    #code below sets up the plots
    ax1.clear()
    ax1.set_xlabel('k vector (as fraction of G)');
    ax1.set_ylabel('Electron Energy');
    ax1.set_title('Band Structure')
    
    for i in range(0,n_G+1):
        ax1.plot(k_values[:,i],Energy_Eigenvalues[:,i],'b')
        ax1.plot(k_values[:,i]+1,Energy_Eigenvalues[:,i],'g')
        ax1.plot(k_values[:,i]-1,Energy_Eigenvalues[:,i],'r')
   
    ax2.clear()
    ax2.set_xlabel('Distance x');
    ax2.set_ylabel('Potential U(x)');
    ax2.set_title('Shape of Potential U(x)')
    ax2.set_yticks(ticks=[-1,0,1]);
    ax2.plot(x[int((n/2-1/2)*fs)-2:int((n/2+1/2)*fs)+2], V[int((n/2-1/2)*fs)-2:int((n/2+1/2)*fs)+2])
    
    ax3.clear()
    ax3.plot(k, Vq.real)
    ax3.set_xlim(0,10)
    ax3.set_xticks(ticks=np.arange(0, 10, step=1));
    ax3.set_yticks(ticks=[0]);
    ax3.set_title('Fourier transform of Potential U(q)')
    ax3.set_xlabel('q vector (as fraction of G)');
    ax3.set_ylabel('Relative Intensity of Uq');
    ax3.grid
    
    x = np.linspace(-1.5*n_k,1.5*n_k,100)
    plz=t*(x/(n_k))**2 + np.ones(x.shape[0])*Energy_Eigenvalues[int(n_k/2+1)][0]
    x = x/n_k
    ax1.plot(x,plz,'y')
    ax1.set_ylim(Energy_Eigenvalues[int(n_k/2+1)][0]-0.1,(Energy_Eigenvalues[int(n_k/2+1)][7]+Energy_Eigenvalues[int(n_k/2+1)][8])/2)

    plt.show()
In [3]:
#code below is mostly identical to the one above, but with fourier transforms of potentials removed
#changing values of t and u doesn't change the general shape of potential, hence we only need to update the band structure
def interactive_bands_tu(t,u):
    '''
    :param k: wave number
    '''
    global t_g, u_g, U_type_g, n_G, n_k
    global ax1, ax2, ax3
    t_g = t #this time round we update u and t hence we pass these to global variables
    u_g = u
    U_type = U_type_g #read the current potentail type
    mat=np.zeros((n_G+1,n_G+1)) #predefine the size of the Hamiltonian matrix
    k_values=np.zeros((n_k+1,n_G+1)) #predefine the size of variables which will be used for band plots
    Energy_Eigenvalues=np.zeros((n_k+1,n_G+1))
    
   
    
    if (U_type == 'Sine Wave'):
        for i in range(n_k+1):
            mat=np.zeros((n_G+1,n_G+1))
            U=np.ones(n_G+1)*(-u) #define fourier components of the potential. As potential is cos() hence only G'-G = +-1 compontents will be present
            K=np.ones(n_G+1) #define dimensions of kinetic energy vector
            for j in range(n_G+1):
                K[j] = t*(i/n_k - 1/2 + j - n_G/2)**2 #for every
            mat=mat+np.diag(U,-1)[:-1,:-1]+np.diag(U,1)[:-1,:-1]+np.diag(K,0) 
            #put the kinetic terms on the diagonal of the Hamiltonian, and add the potential as off diagonal elements
            Energy_Eigenvalues[i]=np.linalg.eigvalsh(mat) #energy levels will be the eigenvalues of hamiltonian
            k_values[i]=np.ones(n_G+1)*(i-n_k/2)/n_k #keep track of which k value is the eigenvalues for

    elif (U_type == 'Series of Delta Functions'):
        for i in range(n_k+1):
            K=np.ones(n_G+1)
            for j in range(n_G+1):
                K[j] = t*(i/n_k - 1/2 + j - n_G/2)**2
            mat=np.ones((n_G+1,n_G+1))*(u) #in this case, all G'-G values are u, hence we just define all matrix elements to be u
            mat=mat+np.diag(K,0) #add kinetic term to the diagonal
            Energy_Eigenvalues[i]=np.linalg.eigvalsh(mat) #energy levels will be the eigenvalues of hamiltonian
            k_values[i]=np.ones(n_G+1)*(i-n_k/2)/n_k 
            
    elif (U_type == 'Square Wave'):
        mat=np.zeros((n_G+1,n_G+1))
        for i in range(n_G+1):
            for j in range(n_G+1):
                if ((i - j) %2!=0):
                    mat[i][j] = (-1)**(int((i-j+1)/2))*u/(i-j) 
        #in this case, the fourier coefficients will be non-zero only if G'-G is odd, are proportional to 1/(G'-G)
        #and vary in sign, -1 if (G'-G)/2pi = = 1 (mod 4), 1 if (G'-G)/2pi = -1 (mod 4) [a=1 for convinience]
        for i in range(n_k+1):
            mat_sq = np.copy(mat)
            K=np.ones(n_G+1)
            for j in range(n_G+1):
                K[j] = t*(i/n_k - 1/2 + j - n_G/2)**2
            mat_sq=mat_sq+np.diag(K,0) #add kinetic energy to the hamiltonian
            Energy_Eigenvalues[i]=np.linalg.eigvalsh(mat_sq) #energy levels will be the eigenvalues of hamiltonian
            k_values[i]=np.ones(n_G+1)*(i-n_k/2)/n_k 
    

    ax1.clear()
    ax1.set_xlabel('k vector (as fraction of G)');
    ax1.set_ylabel('Electron Energy');
    ax1.set_title('Band Structure')
    
    for i in range(0,n_G+1):
        ax1.plot(k_values[:,i],Energy_Eigenvalues[:,i],'b')
        ax1.plot(k_values[:,i]+1,Energy_Eigenvalues[:,i],'g')
        ax1.plot(k_values[:,i]-1,Energy_Eigenvalues[:,i],'r')

    x = np.linspace(-1.5*n_k,1.5*n_k,100)
    plz=t*(x/(n_k))**2 + np.ones(x.shape[0])*Energy_Eigenvalues[int(n_k/2+1)][0]
    x = x/n_k
    ax1.plot(x,plz,'y')
    ax1.set_ylim(Energy_Eigenvalues[int(n_k/2+1)][0]-0.1,(Energy_Eigenvalues[int(n_k/2+1)][7]+Energy_Eigenvalues[int(n_k/2+1)][8])/2)

An interactive look¶

In the following block, there is a simple example on how band structure emerges in lattices. The potential function $U$ is for a lattice with atomic positions $x_i \in \mathbb{Z}$. $U$ takes one of three forms:

  1. Sine wave: $U(x)=-2u*cos(2\pi x)$,

  2. Series of Delta functions centred on atomic positions: $U(x)=\sum_{x_i} u*\delta(x - x_i)$

  3. Square wave: $U(x)= -u$ for $x_i - \frac{1}{2} < x \le x_i + \frac{1}{2}$ and $U(x)= u$ otherwise

By tuning the parameter u, one can travel between the periodic free electronic limit and the discrete energy levels. u is defined in such a way that at u=1, it's contribution to the Hamiltonian will be similar to the kinetic energy contribution. One can see that, at the limit $1<<u$, the spectrum resembles a series of discrete energy levels, which are rarely dependent on the wave vector $k$. Such a spectrum reflect the fact that the eigenstates are mostly localized by the strong potential. On the other hand, if we take $u<<1$, the spectrum then looks like a series of quadratic functions. This is consistent with the fact that the system is periodic and the electrons are nearly free.

In between the two cases, there then arrives our band structures: the spectrum forms several energy bands while each band has some dispersion.

In [4]:
plt.clf
figs = plt.figure();
ax1 = plt.subplot(1,9,(1,3));

ax2 = plt.subplot(3,9,(6,9));

ax3 = plt.subplot(2,9,(15,18));


interact(interactive_bands_tu,t=fixed(1),u=widgets.FloatSlider(min=0.00, max=1.00, step=0.05, value = u_g));
# input t=widgets.FloatSlider(min=0.05, max=1.00, step=0.05, value = t_g) instead of t=fixed(1) to have control over KE contribution
interact(interactive_bands_V, U_type=widgets.RadioButtons(
    value= U_type_g, 
    options=['Sine Wave', 'Series of Delta Functions', 'Square Wave'], 
    description='Potential:'
));

interactive(children=(FloatSlider(value=0.5, description='u', max=1.0, step=0.05), Output()), _dom_classes=('w…
interactive(children=(RadioButtons(description='Potential:', options=('Sine Wave', 'Series of Delta Functions'…
In [ ]: