I am trying to figure out the numerical interpretation of density of states for a fermionic system under a periodic potential.
The equation for the density of states reads
$$
DOS(E) = \sum_{k \in BZ, n} \delta(E-E_n(k)),
$$
where $E_n(k)$ are the eigenvalues of the particular Hamiltonian matrix I am solving. I would like to use the Cauchy/Lorentzian approximation of the Delta function such that the first equation now becomes
$$
DOS(E) = \frac{1}{\pi} \lim_{\epsilon\to 0} \sum_{k \in BZ, n} \frac{\epsilon}{(E-E_n(k))^2+\epsilon^2}.
$$
From here on, I am confused on how to interpret the second equation numerically. I have the respective eigenvalues of the Hamiltonian, but I don't know how to obtain the DOS using $E$. How do I include $E$ in my code? Discretizing E to me means that I grab a certain energy window around a certain value $E$, but I don't know how to structure it, or if needs to be an array, a grid…or something else. If E should be a grid, should it be a grid between the minimum and maximum values of the energy eigenvalues?
EDIT:
Hey all. After pondering about Murali's answer I have come up with a pseudo code that is rather bad but I would like to know if I am going in the right direction.
I basically coded a function for the Lorentzian broadened delta function like so:
def delta_l(x):
return (1/np.pi)*(epsilon/(epsilon**2 + x**2))
def dos(Egrid,Eigen):
DOS = np.zeros((AllK,1))
for j in range(Allk):
DOS[j] = (1/AllK)*sum([delta(Egrid[j]-Eigen[i]) for i in range(np.shape(Eigen)[0])])
return DOS
Here epsilon is given a value of 0.1 just to test. The eigenvectors of the Hamiltonian were obtained by inputting points of the FBZ:
AllK = len(np.arange(0, 1, 0.01)) * len(np.arange(0, 1, 0.01))
E = np.zeros((AllK,4*n), float)
count = 0
for m in np.arange(0, 1, 0.01):
for f in np.arange(0, 1, 0.01):
kx = (m-f) * np.sqrt(3)/2
ky = (m+f) * 3.0/2 - 1
E[count] = Hamiltonian(kx*Kmag, ky*Kmag)
count = count + 1
import pandas as pd
EinBZ = E.flatten()
So then I get all the eigenvalues of the FBZ in this array. Am I going in the right direction?
Best Answer
I think you did not correctly understand what density of states (DOS) mean. DOS is a probability density function (PDF). As Andrew pointed out, it takes energy as input and returns the number of states for a given energy.
You cannot discretize $E$ as they are not eigenvalues of any observable. It is the input parameter, and discretizing it simply does not make any sense. The values of $E_n(k)$ are discretized as they are eigenvalues of the electronic hamiltonian.
If you consider the 1st equation in your question,
$$ DOS(E) = \sum_{k \in BZ, n} \delta(E-E_n(k)), $$ For Energies $E\neq E_n(k), DOS(E) = 0$. However, this is not what you observe in experiments. Typically we see some finite DOS for energies not equal to eigenstates due to Heisenberg uncertainty. To account for this, we add a small finite electronic broadening parameter ($\epsilon$) as shown in the second equation of your question.
To compute the DOS, you take E as a parameter which can take any value and fix $\epsilon$, then compute the double summation over Brillouin zone (BZ) and bands (n) which are the eigenvalues of your hamiltonian. Summing over BZ is simply summation over k points in the Brillouin zone and divide the obtained sum by the total number of k points. Choose a reasonable k point grid and make sure it is converged. Have a look at the following link (http://www.iiserpune.ac.in/~smr2626/hands_on/week1/july1/bzsums_mastani.pdf) if you don't have any idea about BZ summation
Pseudo code: