Quantum Band Structure Modeling- From 1D chains to Graphene

Author: Aneesh Jonelagadda          Date: 04/2023 (Edited 2025)

Band diagrams represent the energy levels and behavior of electrons at various symmetry points in an atomic lattice. Calculating the energy states and band structure requires solving the Schrödinger equation in quantum mechanics, which is not tractable for multi-atom and multi-electron systems such as lattices or physical materials. Condensed-matter physics introduces several levels of approximating the Schrödinger Equation such that it can be solved analytically and/or numerically, and one of these approximations is called the "tight-binding" model. The tight binding model describes an electron by treating its wave function as a combination of simpler localized atomic orbital wavefunctions, and uses hopping integrals to represent the probability of an electron moving between atomic sites.

I will use the tight binding model in this project to model several condensed matter systems. This will include:

1. Face-Centered-Cubic (first with nearest neighbors, and then including next-nearest neighbors)
2. Body-Centered-Cubic (up to next-nearest neighbors)
3. Diatomic Linear Chain of H-Li
4. Graphene

1 and 2 will be modeled fully numerically by lattice vector-based expansion. 3 and 4 will be first modeled analytically then plotted numerically. Additional commentary will be provided to explain steps and provide context for how the calculations and modeling relate back to physical properties, and commentary will be geared towards a graduate audience familiar with quantum physics.

Initial Step 1: How do we plot a band diagram?

We note that we are plotting against a panel of directions in a lattice, and thus create a function to create an panel with helper boundary points that give us the plotting indices of the high symmetry points so that we can mark these on the final band diagram.

The main objective output of the panel function is a discretized k-space grid across the input panel directions with spacing corresponding to the physically realistic distances between the high-symmetry points. We can use this grid to probe the energy in a nice pythonic vectorized manner (plus it's useful for plotting).

import numpy as np
from IPython.display import Latex
import matplotlib.pyplot as plt
plt.rcParams['axes.xmargin'] = 0
plt.rcParams['axes.ymargin'] = 0

#Defines helper function which returns the squared euclidean distance between two points
def dist(a,b):
    #If 3d vector:
    if len(a)==3:
        return(((a[0]-b[0])**2+(a[1]-b[1])**2+(a[2]-b[2])**2)**0.5)
    #If 2d vector
    elif len(a)==2:
        return(((a[0]-b[0])**2+(a[1]-b[1])**2)**0.5)

#Defines the panel function, creates the backbone of the band diagrams
def panel(point_sequence,a=1,grid_density = 250): #a default is 1 Angstrom
    '''
    Input: k-space point sequence of the panel
    Outputs: 
    k_grid: Parametrized k-grid across the panel
    k_grid_xaxis: x axis in k-space of the panel
    boundary_indices: Locations in index along the parametrized grid which correspond to the input k-space points
    '''
    pi=np.pi
    
    #Define the grid-point index to k space conversion factor. 
    #Note we define the grid_density to be # of grid points per (2pi/a) unit length
    scale_factor = grid_density/(2*pi/a)
    
    k_grid = [] #Initializes our empty grid
    boundary_indices = [0]
    boundary_index = 0
    #Forms the direction pairs
    for i in range(0,len(point_sequence)-1):
        kspace_dir_point1 = point_sequence[i]
        kspace_dir_point2 = point_sequence[i+1] #Now we have the points for which to interpolate
        
        #Calculates the distance between the two points, we will need this for our plot spacings
        kspace_dir_dist = dist(kspace_dir_point1,kspace_dir_point2)
        
        #Gets the relative length of the direction compared to the gamma-X direction, uses base grid density
        rel_length = int((kspace_dir_dist*scale_factor))
        boundary_index += rel_length #Adds the relative length to boundary index to keep track of where we are!
        
        #Generates k-space subgrid corresponding to the direction with scaled length
        kspace_subgrid = np.linspace(kspace_dir_point1,kspace_dir_point2,rel_length)
        
        #Adds subgrid to the grid, adds boundary
        k_grid.extend(kspace_subgrid)
        boundary_indices.append(boundary_index)
        
    #Converts the x axis of our panel to k-space points instead of grid points.
    panel_xaxis = np.arange(0,len(k_grid))/scale_factor
    return(np.array(k_grid),panel_xaxis,boundary_indices)

Initial Step 2: Energy function

We now create the meat of this, which is the generic energy function. For now this will correspond to the tight-binding energy FOR 1S ORBITALS for Bloch eigenfunctions:

$$E_\vec k = E(\vec k) = E_{n=1s} - \beta_s - \sum_{\vec R}t_{\vec R}e^{i\vec k \cdot \vec R}$$

For the sake of systems #1 and #2, we assume $E_{n=1s} - \beta_s = 0$ and $t_{\vec R}$ the nearest-neighbor surface hopping integral = 1 eV

The above equation can be changed to, for nearest neighbors in 3D: $$E_\vec k = E_{n=1s} - \beta_s - \frac{1}{2} t_{\vec R} \sum_{\vec R} \cos(\vec k \cdot \vec R)$$

def tbind_energy_1s_nn(k_vec, R_vectors, E_s = 0, beta_s = 0, t_R = 1):
    '''
    k_vec: k vector (from grid)
    R_vectors: nearest neighbor lattice vectors
    t_R = nearest-neighbor surface hopping integral
    '''
    energy = 0
    #Sums across each lattice vector
    for R_vec in R_vectors:
        energy += E_s - beta_s - 0.5*t_R*np.cos(np.dot(k_vec,R_vec))
    #Returns the final energy
    return(energy)

System 1 Modeling: Face-Centered-Cubic (fcc)

Face-centered cubic lattice structure is a common crystal structure seen in Copper, Aluminum, and Gold for example. Here's a nice visualization, using Crystalmaker

k-space point defining

We define the high symmetry points in k-space for fcc mathematically, and set a = 4 Å. We can of course change this whenever.

pi = 3.1415927
a = 4 #We are setting a to be 4 Å

Γ = np.array([0,0,0])
X = (2*pi/a)*np.array([1,0,0])
L = (pi/a)*np.array([1,1,1])
K = (2*pi/a)*np.array([0.75,0.75,0])
W = (2*pi/a)*np.array([1,0.5,0])

Now we create our iteration function and run the computation

We first initialize the primary nearest neighbor lattice vectors for fcc. We can see these lattice vectors to the 12 first nearest neighbors visualized below from the central atom:

#Initializes the 'basis vectors' (we adjust the length of each basis to be a/2 later)
nn_vec1 = np.array([1,1,0])
nn_vec2 = np.array([1,0,1])
nn_vec3 = np.array([0,1,1])

#Now we include symmetry and create our 12 nearest neighbor vectors
nn_vecs = []
for i in (-1,1):
    for j in (-1,1):
        for k in (-1,1):
            nn_vecs.append(np.multiply(nn_vec1,np.array([i,j,k])))
            nn_vecs.append(np.multiply(nn_vec2,np.array([i,j,k])))
            nn_vecs.append(np.multiply(nn_vec3,np.array([i,j,k])))
    
nn_vecs = np.array(nn_vecs)*(a/2)

#Now we run our calculation, first by getting our k-grid using the panel function.

#Creates panel based on our point sequence:
grid_density = 250 #Sets number of grid points per 2pi/a unit length in k space
scale_factor = grid_density/(2*pi/a) #Gets a conversion factor from grid index to k-space coords
point_seq = [L,Γ,X,K,Γ,W]
k_grid, panel_xaxis, bounds = panel(point_seq,a=a,grid_density=grid_density)

#Iterates over our k space grid and calculates tight binding energies
energy_band = []
for k in k_grid:
    energy = tbind_energy_1s_nn(k_vec=k, R_vectors=nn_vecs)
    energy_band.append(energy)

#Plots energy band
plt.figure(figsize=[10,7])
plt.plot(panel_xaxis,energy_band,color='black')
plt.ylabel("Energy (eV)",fontsize=13)
plt.title("Tight Binding Band Diagram for FCC (1s orbitals, Nearest-Neighbors only)",fontsize=15)

#Plots the high symmetry points as dotted lines
#Note we need to convert these points from grid coords to k-space by dividing by grid density
for sympoint_index in bounds:
    plt.axvline(sympoint_index/scale_factor,color='black',linestyle='dotted')

#Marks the high symmetry points on x axis, labels them with x ticks
plt.xticks(np.array(bounds)/scale_factor,['L','Γ','X','K','Γ','W'],fontsize=17)
plt.show()

Band-width

The bandwidth is $E_{max} - E_{min}$ = 4 eV - (-12 eV) = 16 eV.

We use the bandwidth to measure roughly how localized the electrons are; if the bandwidth is narrow then electrons are more localized. In this case this would correspond to pretty delocalized electrons.

Effective Mass of Electrons

The effective mass of an electron is highly important to calculate carrier mobility, which has immense impact in conductivity properties.

With a= 4 Å, we use a parabolic approximation around the Γ point.

We can find the derivative $\frac{\partial^2 E}{\partial k^2}$ from the parabolic fit $E = A_0 k^2 + A_1 k + A_2$ , and see this is simply $\frac{\partial^2 E}{\partial k^2} = 2A_0$. Thus, given definition of effective mass $$m^* = \hbar ^2/\left(\frac{\partial^2 E}{\partial k^2}\right)$$ we get: $$m^* = \hbar ^2/(2A_0)$$

#Selects band close to gamma point
e_band_gamma = energy_band[bounds[1]-50:bounds[1]+50]
x_band_gamma = panel_xaxis[bounds[1]-50:bounds[1]+50]
#Does polynomial fit
parabolic_fit = np.poly1d(np.polyfit(x_band_gamma,e_band_gamma,2))
print("Fitted Polynomial:")
print(parabolic_fit)
#Plots the fit of the small portion of the band for a sanity check
plt.xlabel("k (1/Å)")
plt.ylabel("Energy (eV)")
plt.plot(x_band_gamma,e_band_gamma,color='black',label='Calculated Band')
plt.plot(x_band_gamma,parabolic_fit(x_band_gamma),color='red',linestyle='dotted',label='Fitted Parabola')
plt.legend()
plt.show()

Fitted Polynomial:
       2
15.37 x - 41.61 x + 16.17

Thus since $A_0 = 15.37 $ eV $Å^2$, the effective mass is $m^* = \hbar ^2/(2A_0)$ = 2.258e-31 kilograms

FCC with Next-nearest-neighbors

There are now new lattice vectors R' that can take us to areas on the lattice where the overlap integral may be significant. These R' will point to the next-nearest neighbors, and the surface-hopping integral at these sites will be different, called t'. This gives us:

$$E_\vec k = E_{n=1s} - \beta_s - \frac{1}{2} t_{\vec R} \sum_{\vec R} \cos(\vec k \cdot \vec R)- \frac{1}{2} t'_{\vec R} \sum_{\vec R'} \cos(\vec k \cdot \vec R') $$

For this problem we assume $t_{\vec R}$ = 1 eV and the next-nearest-neighbor surface hopping integral is $t'_{\vec R}$ = 0.5 eV. We thus modify our energy function:

def tbind_energy_1s_nnn(k_vec, R_vectors, R_prime_vectors, E_s = 0, beta_s = 0, t_R = 1,t_prime_R = 0.5):
    '''
    k_vec: k vector (from grid)
    R_vectors: nearest neighbor lattice vectors
    R_prime_vectors: next-nearest neighbor lattice vectors
    t_R = nearest-neighbor surface hopping integral
    t_prime_R = next-nearest-neighbor surface hopping integral
    '''
    energy = 0
    #Sums across each nearest neighbor lattice vector
    for R_vec in R_vectors:
        energy += E_s - beta_s - 0.5*t_R*np.cos(np.dot(k_vec,R_vec))
    #Sums across each next-nearest neighbor lattice vector
    for R_prime_vec in R_prime_vectors:
        energy += -0.5*t_prime_R*np.cos(np.dot(k_vec,R_prime_vec))
    #Returns the final energy
    return(energy)

Now we define the next-nearest lattice vectors. This is easy; these are just the 'simple cubic' directions (along each axis) and there are 6 of these. After we do this, we proceed with the same calculation as before.

#Initializes the 'next-nearest-neighbor basis vectors'
nnn_vec1 = np.array([1,0,0])
nnn_vec2 = np.array([0,1,0])
nnn_vec3 = np.array([0,0,1])

#Now we include symmetry and create our next-nearest neighbor vectors
nnn_vecs = []
for i in (-1,1):
    for j in (-1,1):
        for k in (-1,1):
            nnn_vecs.append(np.multiply(nnn_vec1,np.array([i,j,k])))
            nnn_vecs.append(np.multiply(nnn_vec2,np.array([i,j,k])))
            nnn_vecs.append(np.multiply(nnn_vec3,np.array([i,j,k])))
    
nnn_vecs = np.array(nnn_vecs)*(a)

#Initializes the 'nearest neighbor basis vectors' (This is the same as before)
nn_vec1 = np.array([1,1,0])
nn_vec2 = np.array([1,0,1])
nn_vec3 = np.array([0,1,1])

#Now we include symmetry and create our nearest neighbor vectors
nn_vecs = []
for i in (-1,1):
    for j in (-1,1):
        for k in (-1,1):
            nn_vecs.append(np.multiply(nn_vec1,np.array([i,j,k])))
            nn_vecs.append(np.multiply(nn_vec2,np.array([i,j,k])))
            nn_vecs.append(np.multiply(nn_vec3,np.array([i,j,k])))
    
nn_vecs = np.array(nn_vecs)*(a/2)

#Now we run our calculation, first by getting our k-grid using the panel function.

#Creates panel based on our point sequence:
grid_density = 250 #Sets number of grid points per unit length in k space
scale_factor = grid_density/(2*pi/a) #Gets a conversion factor from grid index to k-space coords
point_seq = [L,Γ,X,K,Γ,W]
k_grid, panel_xaxis, bounds = panel(point_seq,a=a,grid_density=grid_density)

#Iterates over our k space grid and calculates tight binding energies
energy_band = []
for k in k_grid:
    energy = tbind_energy_1s_nnn(k_vec=k, R_vectors=nn_vecs, R_prime_vectors=nnn_vecs)
    energy_band.append(energy)

#Plots energy band
plt.figure(figsize=[10,7])
plt.plot(panel_xaxis,energy_band,color='black')
plt.ylabel("Energy (eV)",fontsize=13)
plt.title("Tight Binding Band Diagram for FCC (1s orbitals, Next Nearest-Neighbors Included)",fontsize=15)

#Plots the high symmetry points as dotted lines
#Note we need to convert these points from grid coords to k-space by dividing by grid density
for sympoint_index in bounds:
    plt.axvline(sympoint_index/scale_factor,color='black',linestyle='dotted')

#Marks the high symmetry points on x axis, labels them with x ticks
plt.xticks(np.array(bounds)/scale_factor,['L','Γ','X','K','Γ','W'],fontsize=17)
plt.show()

System 2 Modeling: Body-Centered Cubic (bcc)

Here's a nice visualization of a body-centered cubic lattice:

To model this, we first define new symmetry points for our panel.

pi = 3.1415927
a = 4 #We are setting a to be 4 angstrom

Γ = np.array([0,0,0])
N = (2*pi/a)*np.array([0.5,0.5,0])
P = (2*pi/a)*np.array([0.5,0.5,0.5])
H = (2*pi/a)*np.array([0,1,0])

Now we define the new nearest-neighbor and next-nearest-neighbor vectors, and proceed with the problem. Note we plot both the NN band and NNN bands in same panel. For BCC our nearest-neighbor vectors are [+/- a/2,+/- a/2,+/- a/2]. The next-nearest-neighbors are the same as in FCC.

#Initializes the 'next-nearest-neighbor basis vectors'. Luckily for bcc this is same as FCC.
nnn_vec1 = np.array([1,0,0])
nnn_vec2 = np.array([0,1,0])
nnn_vec3 = np.array([0,0,1])

#Now we include symmetry and create our next-nearest neighbor vectors
nnn_vecs = []
for i in (-1,1):
    for j in (-1,1):
        for k in (-1,1):
            nnn_vecs.append(np.multiply(nnn_vec1,np.array([i,j,k])))
            nnn_vecs.append(np.multiply(nnn_vec2,np.array([i,j,k])))
            nnn_vecs.append(np.multiply(nnn_vec3,np.array([i,j,k])))
    
nnn_vecs = np.array(nnn_vecs)*(a)

#Initializes the 'nearest neighbor basis vectors'. For BCC this is just [+/- a/2,+/- a/2,+/- a/2] with 
nn_vec = np.array([1,1,1])

#Now we include symmetry and create our nearest neighbor vectors
nn_vecs = []
for i in (-1,1):
    for j in (-1,1):
        for k in (-1,1):
            nn_vecs.append(np.multiply(nn_vec,np.array([i,j,k])))
    
nn_vecs = np.array(nn_vecs)*(a/2)

#Now we run our calculation, first by getting our k-grid using the panel function.

#Creates panel based on our point sequence:
grid_density = 250 #Sets number of grid points per unit length in k space
scale_factor = grid_density/(2*pi/a) #Gets a conversion factor from grid index to k-space coords
point_seq = [Γ,N,P,H]
k_grid, panel_xaxis, bounds = panel(point_seq,a=a,grid_density=grid_density)

#Iterates over our k space grid and calculates tight binding energies
energy_band_nn = []
energy_band_nnn = []
for k in k_grid:
    energy_nn = tbind_energy_1s_nn(k_vec=k, R_vectors=nn_vecs)
    energy_nnn = tbind_energy_1s_nnn(k_vec=k, R_vectors=nn_vecs, R_prime_vectors=nnn_vecs)
    energy_band_nn.append(energy_nn)
    energy_band_nnn.append(energy_nnn)
    
#Plots energy band
plt.figure(figsize=[10,7])
plt.plot(panel_xaxis,energy_band_nn,color='black',label='Nearest-Neighbor')
plt.plot(panel_xaxis,energy_band_nnn,color='green',label='Next-Nearest-Neighbor')
plt.ylabel("Energy (eV)",fontsize=13)
plt.title("Tight Binding Band Diagram for BCC (1s orbitals, Next Nearest-Neighbors Included)",fontsize=15)

#Plots the high symmetry points as dotted lines
#Note we need to convert these points from grid coords to k-space by dividing by grid density
for sympoint_index in bounds:
    plt.axvline(sympoint_index/scale_factor,color='black',linestyle='dotted')

#Marks the high symmetry points on x axis, labels them with x ticks
plt.xticks(np.array(bounds)/scale_factor,['Γ','N','P','H'],fontsize=17)
plt.legend()
plt.show()

System 3 Modeling: H-Li Diatomic Linear Chain

We now model a diatomic linear chain of alternating Hydrogen and Lithium atoms.

In problems 1 and 2 we solved the system fully numerically, meaning we relied on computational methods to numerically expand the $\sum_{\vec R} \cos(\vec k \cdot \vec R)$ and $\sum_{\vec R'} \cos(\vec k \cdot \vec R')$ across the possible k-space grid $\vec k$ and lattice vectors $\vec R$ in $$E_\vec k = E_{n=1s} - \beta_s - \frac{1}{2} t_{\vec R} \sum_{\vec R} \cos(\vec k \cdot \vec R)- \frac{1}{2} t'_{\vec R} \sum_{\vec R'} \cos(\vec k \cdot \vec R')$$

This approach gives concrete results and can be definitely applied to both this system and Graphene in System 4. However, for problems 3 and 4 we choose to derive the dispersion analytically and derive closed-form expressions by directly expanding the tight-binding Hamiltonian. This produces explicit formulas that reveal how the energies depend on parameters and symmetry. The analytical route provides more physical insight while the numerical route lays a foundation which can be scaled to any lattice in theory by changing the lattice basis vectors and k-space grid.

We begin from the time-independent Schrödinger equation:

$$ \hat{H}\,\Psi_{\mathbf{k}} = E_{\mathbf{k}}\,\Psi_{\mathbf{k}} . $$

For a diatomic linear chain, the basis contains two types of atoms per unit cell, so the wavefunction coefficients are written as a two-component vector $(c_A, c_B)$. The solution for matrix form of the system of Schrödinger equations thus becomes (for non-trivial solutions):

$$ \det \begin{bmatrix} H_{AA} - E & H_{AB} \\ H_{BA} & H_{BB} - E \end{bmatrix} = 0 $$

Next, we substitute the Hamiltonians using the tight-binding approximation with decoupling between on-site energies and site-hopping integrals with phase factors:

$$ H_{AA} = H_{BB} = \epsilon \quad\text{(on-site energies)}, $$

$$ H_{AB} = H_{BA}^* = -t\sum_m e^{ik\cdot R_m} \quad\text{(site-hopping with phase factors)}. $$

That substitution turns the determinant into

$$ (\epsilon - E)^2 - t^2 \left|\,1 + e^{ik\cdot a_1} + e^{ik\cdot a_2}\,\right|^2 = 0 , $$

and when expanded, this produces the cosine terms that depend explicitly on $k_x, k_y$.
Those cosines reflect the lattice symmetry, and solving for $E$ gives the analytical dispersion relation.

From the pen+paper algebra, we get: $$ E(k) = \frac{1}{2}(E_A + E_B) + 2t_2\cos(ka) \pm \sqrt{\frac{1}{4}(E_A-E_B)^2+ 4t_1^2\cos^2(ka/2)}$$

For the first plot we assume $t_1$ = 0.5 eV, $t_2$ = 0.1 eV, $E_A = E_B = 0$, and a= 10 Å. For the second plot we assume $E_A = 0.8 eV$ and $E_B = 0.6 eV$, from the values of H and Li found in https://doi.org/10.48550/arXiv.1112.4638

a = 10 #Assume a= 10 Å for the sake of plotting
E_A,E_B = 0,0
t_1,t_2 = 0.5,0.1
#Generates grid according to 1st brillouin zone
k_grid = np.linspace(-pi/a,pi/a,250)

#Lower band energy
def energy_dc_lower(k,E_A,E_B,t_1,t_2,a):
    return(0.5*(E_A+E_B)+2*t_2*np.cos(k*a)-np.sqrt(0.25*(E_A-E_B)**2+4*t_1**2*(np.cos(k*a*0.5))**2))
#Upper band energy
def energy_dc_upper(k,E_A,E_B,t_1,t_2,a):
    return(0.5*(E_A+E_B)+2*t_2*np.cos(k*a)+np.sqrt(0.25*(E_A-E_B)**2+4*t_1**2*(np.cos(k*a*0.5))**2))

#Plots both bands
plt.figure(figsize=[8,6])
plt.plot(k_grid,energy_dc_lower(k_grid,E_A=0,E_B=0,t_1=t_1,t_2=t_2,a=a),color='black')
plt.plot(k_grid,energy_dc_upper(k_grid,E_A=0,E_B=0,t_1=t_1,t_2=t_2,a=a),color='black')
plt.xlabel("k (1/Å)")
plt.ylabel("Energy (eV)")
plt.title("Tight Binding Bands for Diatomic Chain, E_A=E_B=0")
plt.show()

#Plots both bands for if E_A and E_B are different and nonzero.
#We choose E_A to be 0.8 eV and E_B to be 0.6 eV.
plt.figure(figsize=[8,6])
plt.plot(k_grid,energy_dc_lower(k_grid,E_A=0.8,E_B=0.6,t_1=t_1,t_2=t_2,a=a),color='black')
plt.plot(k_grid,energy_dc_upper(k_grid,E_A=0.8,E_B=0.6,t_1=t_1,t_2=t_2,a=a),color='black')
plt.xlabel("k (1/Å)")
plt.ylabel("Energy (eV)")
plt.title("Tight Binding Bands for Diatomic Chain, E_A=0.8 eV E_B=0.6 eV")
plt.show()

Effective Mass

For this, we simply follow the same procedure as #1 regarding the parabolic approximation. Recall the relation: $$m^* = \hbar ^2/(2A_0)$$ where $A_0$ is our 2nd-order coefficient of parabolic fit.

#Selects band close to center
k_grid_small = k_grid[75:175]
e_dc_lower_small = energy_dc_lower(k_grid_small,E_A=0,E_B=0,t_1=t_1,t_2=t_2,a=a)
#Does polynomial fit
parabolic_fit = np.poly1d(np.polyfit(k_grid[75:175],e_dc_lower_small,2))
print("Fitted Polynomial:")
print(parabolic_fit)
#Plots the fit of the small portion of the band for a sanity check
plt.xlabel("k (1/Å)")
plt.ylabel("Energy (eV)")
plt.plot(k_grid_small,e_dc_lower_small,color='black',label='Calculated Band')
plt.plot(k_grid_small,parabolic_fit(k_grid_small),color='red',linestyle='dotted',label='Fitted Parabola')
plt.legend()
plt.show()

Fitted Polynomial:
       2
3.236 x - 2.439e-15 x - 0.8011

Thus since $A_0 = 3.236 $ eV $Å^2$, the effective mass is $m^* = \hbar ^2/(2A_0)$ = 1.07e-30 kilograms

Since there is a pretty visually obvious band-gap at all points in the 1st BZ, the material will be an insulator since there is no overlap in the bands.

System 4 Modeling: Graphene

Graphene is an interesting 2D material made up of tiled Carbon hexagons. It exhibits exotic electron and phonon behavior due to various quirks, a few of which can be explained by examining the band diagram.

We proceed similarly to System 3, but first we calculate the reciprocal lattice vectors for Graphene and define the symmetry. Then we define the system of tight-binding Hamiltonians, and expand. We will use Dirac bra-ket notation to directly represent the wavefunctions and site-hopping e.g $\langle \phi_{\text{2}_{p_{z_2}}} | \hat H |\Psi_k (r) \rangle $ and $\langle \phi_{\text{1}} | \hat H |\phi_{\text{1}}\rangle $

We get from our paper and pen calculations the dispersion relation for graphene in our tight-binding approximation model. While we could do this purely computationally in a 'lattice-agnostic' way using the lattice vectors and the code handling the expansion of the sum across NNs implicitly, we step through analytically in the paper+pen calculations to get a nice analytical form of the dispersion relation:

$$E(\vec k) = \epsilon \pm t\sqrt{1+4\cos(\frac{\sqrt 3 k_x a}{2})\cos(\frac{k_y a}{2}) + 4\cos^2(\frac{k_y a}{2})} $$

Now we define an energy function in python for this: (we assume $\epsilon = 0$ eV, t=3 eV, a = 2.46 Å)

def energy_graphene(k,a=2.46,epsilon=0,t=3):
    #Splits the minus and plus into energy lower and energy upper, returns both
    energy_lower = epsilon-t*np.sqrt(1+4*np.cos(np.sqrt(3)*k[0]*a/2)*np.cos(k[1]*a/2)+4*(np.cos(k[1]*a/2))**2)
    energy_upper = epsilon+t*np.sqrt(1+4*np.cos(np.sqrt(3)*k[0]*a/2)*np.cos(k[1]*a/2)+4*(np.cos(k[1]*a/2))**2)
    return(energy_lower,energy_upper)

Graphene high symmetry points and panel

We need to re-define some important high-symmetry points for our panel.

pi = 3.1415927
a = 2.46 #Å

Γ = np.array([0,0])
M = (2*pi/a)*np.array([1/np.sqrt(3) , 0])
K = (2*pi/a)*np.array([1/np.sqrt(3), (-1/3)])

Now we create the panel and plot.

#Creates the panel, first defines point sequence
point_seq = [Γ,M,K,Γ]

#Defines grid density and scale factor (grid dens is # grid points per 2pi/a unit length in k space)
grid_density = 250
scale_factor = grid_density/(2*pi/a)

#Creates panel, evaluates energy at each k-space grid point
k_grid, panel_xaxis, bounds = panel(point_seq,a=a,grid_density=grid_density)
band_lower = []
band_upper = []
for k in k_grid:
    e_lower,e_upper = energy_graphene(k)
    band_lower.append(e_lower)
    band_upper.append(e_upper)

#Plots
plt.figure(figsize=[10,7])
plt.plot(panel_xaxis,band_lower,color='blue',label='Lower Band')
plt.plot(panel_xaxis,band_upper,color='black',label='Upper Band')
plt.ylabel("Energy (eV)",fontsize=13)
plt.title("Tight Binding Band Diagram for Graphene",fontsize=15)

#Plots the high symmetry points as dotted lines
#Note we need to convert these points from grid coords to k-space by dividing by grid density
for sympoint_index in bounds:
    plt.axvline(sympoint_index/scale_factor,color='black',linestyle='dotted')

#Marks the high symmetry points on x axis, labels them with x ticks
plt.xticks(np.array(bounds)/scale_factor,['Γ','M','K','Γ'],fontsize=17)
plt.legend()
plt.show()

This is a really cool and familiar diagram. We can see that there is no band gap at the K point. The bandgap opens up in the K-Γ and K-M directions.

To discuss the fermi level, we note that each Carbon atom in graphene will contribute exactly 1 electron. This makes the lower band completely filled. Thus, the fermi energy $E_f$ will correspond to the highest energy of the lower band, which will be the intersection at the K point.

There is a linear dispersion (conical intersection) around the fermi level and thus around the K point. We can approximate this using a piecewise linear fit. Alternatively we can take a Taylor's series expansion of the dispersion relation around the K point. We use the sum of squared residuals of the linear regression model to show that our model is approximately linear, and a linear fit is justified over a parabolic fit.

#Gets two chunks of the band: one slightly before K, one slightly after K
band_lower_beforeK = band_lower[bounds[2]-25:bounds[2]]
band_lower_afterK = band_lower[bounds[2]:bounds[2]+25]

xaxis_beforeK = panel_xaxis[bounds[2]-25:bounds[2]]
xaxis_afterK = panel_xaxis[bounds[2]:bounds[2]+25]

#Does a linear fit on these two segments
fit_beforeK,resid,rank,sing,rcond = np.polyfit(xaxis_beforeK,band_lower_beforeK,1,full=True)
print("Linear Fit (right BEFORE K, note x=k in inverse angstrom and the equation = E in eV):")
print(np.poly1d(fit_beforeK))
print("Sum of Squared Residuals")
print(np.round(resid[0],4))

fit_afterK,resid,rank,sing,rcond = np.polyfit(xaxis_afterK,band_lower_afterK,1,full=True)
print("\n")
print("Linear Fit (right AFTER K):")
print(np.poly1d(fit_afterK))
print("Sum of Squared Residuals")
print(np.round(resid[0],4))

Linear Fit (right BEFORE K, note x=k in inverse angstrom and the equation = E in eV):
 
5.832 x - 13.49
Sum of Squared Residuals
0.005


Linear Fit (right AFTER K):
 
-6.929 x + 16.09
Sum of Squared Residuals
0.0017

We can see that the sum of squared residuals of a piecewise linear fit is quite low around the K-point. To more robustly show that the fit should be linear, let's examine the coefficients and residuals of a parabolic fit around the same area.

#Gets the before and after chunks together
band_lower_aroundK = band_lower[bounds[2]-25:bounds[2]+25]

xaxis_aroundK = panel_xaxis[bounds[2]-25:bounds[2]+25]

#Does a linear fit on these two segments
parab_fit_aroundK,resid,rank,sing,rcond = np.polyfit(xaxis_aroundK,band_lower_aroundK,2,full=True)
print("Parabolic Fit around K:")
print(np.poly1d(parab_fit_aroundK))
print("Sum of Squared Residuals")
print(np.round(resid[0],4))

Parabolic Fit around K:
        2
-23.39 x + 107.9 x - 124.6
Sum of Squared Residuals
0.7353

We see the sum of squared residuals is 0.7353 for the parabolic fits rather than 0.005 and 0.0017 for the linear fits, which are around 2 orders of magnitude lower! This is pretty good numerical reasoning for the linear behavior of the bands. For even more robustness, let's examine the linear fit right after Γ, where it shouldn't be as linear as near K.

#Gets portion of band right after Γ
band_lower_afterΓ = band_lower[bounds[0]:bounds[0]+25]

xaxis_afterΓ = panel_xaxis[bounds[0]:bounds[0]+25]

fit_afterΓ,resid,rank,sing,rcond = np.polyfit(xaxis_afterΓ,band_lower_afterΓ,1,full=True)
print("Linear Fit (right AFTER Γ):")
print(np.poly1d(fit_afterK))
print("Sum of Squared Residuals")
print(np.round(resid[0],4))

Linear Fit (right AFTER Γ):
 
-6.929 x + 16.09
Sum of Squared Residuals
0.0122

The SSR of the linear is an order of magnitude higher around the Γ point, implying the linear fit is much more appropriate at the K point than other points!

Let's overlay our linear fits around K on top of the band diagram as a nice wrap-up to this numerical analysis:

plt.figure(figsize=[8,6])
plt.title("Fitted Graphene Band-Diagram around K-point")
plt.xlabel("k (1/Å)")
plt.ylabel("Energy (eV)")

#Plots the calculated band
plt.plot(xaxis_aroundK,band_lower[bounds[2]-25:bounds[2]+25],color='black',label='Calculated Band')

#Plots the linear piecewise fits
plt.plot(xaxis_beforeK,np.poly1d(fit_beforeK)(xaxis_beforeK),color='blue',linestyle='dotted',label='Linear Fit')
plt.plot(xaxis_afterK,np.poly1d(fit_afterK)(xaxis_afterK),color='blue',linestyle='dotted')

#Plots the parabolic fits
plt.plot(xaxis_aroundK,np.poly1d(parab_fit_aroundK)(xaxis_aroundK),color='red',
         linestyle='dotted',label='Parabolic Fit')
plt.legend()
plt.show()

NOTE: The flat portion of the calculated band at the very top is an artifact of my k-space discretization into a grid. If I choose a higher grid density, this flat portion would reduce/vanish.

Massless behavior of Electrons in Graphene:

Now that we have demonstrated that the band structure on either side of the K point is approximately linear, we hark back to the effective mass calculation from the dispersion relation:

$$m^* = \hbar ^2/(\frac{\partial^2 E}{\partial k^2})$$

Now we can evaluate this analytically. We are interested in the dispersion around the K-point, and not just on either side, so the earlier linear fits shouldn't be used. Instead we can combine the piecewise linear fits and approximate the dispersion to be an absolute value function centered around $k=k_K$ and E=0 with "slope" $\alpha$ where $k_K$ is the K-point in k-space.

$$E(\vec k) = \alpha |k-k_K|$$

We can express the first derivative in terms of the Heaviside Step Function $\theta$:

$$\frac{\partial E}{\partial k} = \alpha k (\theta(k-k_K) - \theta(-k+k_K))$$

Now we use the fact that the derivative of the Heaviside Step function is a Dirac Delta function and get:

$$\frac{\partial^2 E}{\partial k^2} = 2\alpha \delta(k-k_F)$$

where $\delta$ is the Dirac delta function. We can see that at $k=k_F$, the Dirac delta function will be infinite. Thus the effective mass which is proportional to the inverse of this will be 0 and the excitations, which cross the lower to upper-band at the band intersection which occurs at the K-point, will be 'massless'.

There is a more robust method to prove 0 effective mass near the K point:

Since our E vs k is linear on either side, this will make the effective mass in the first equation diverge, since we are dividing by 0. We can thus use a different formula for the effective mass, and can use the group velocity definition:

$$m^* = \frac{p}{v_g} = \hbar^2k(\frac{\partial E}{\partial k})^{-1}$$

Now we take the approximation that the dispersion relation approximately takes the linear form (using translations to center k=0 around the K point now): $$E(\vec k) \approx \hbar k v_f$$

where $v_f$ is the Fermi velocity of the electrons in graphene. We thus get for the effective mass formula:

$$m^* = \frac{p}{v_g} \approx \frac{p}{v_f} \approx \frac{\hbar k}{v_f}$$

and see when k=0 (we have translated our system to be at K point when k=0), $m^*=0$ and thus the excitations are massless.