# Getting tight binding density of states more accurately

I calculated numerically the density of states (DoS) for the 3-D tightbinding dispersion $\epsilon(k_x,k_y,k_z)=-2t\,(\cos k_x + \cos k_y + \cos k_z)$ and obtained the following plot [$t=1$ has been chosen].

What I did is summing over $k$-points of the lattice Green's function, $$G(k_x,k_y,k_z,\omega)=\frac{1}{\omega-\epsilon(k_x,k_y,k_z)+i0^+}$$ and finding the DoS from its imaginary part : $D(\omega)=-\frac{1}{\pi}\text{Im} \sum_{k_x}\sum_{k_y}\sum_{k_y} G(k_x,k_y,k_z,\omega)$.

One can easily notice that there are noises at low energies. Is there any alternative way to get better result? Like the one shown in a figure from a paper [Ref: arXiv:1207.4014] :

Can there be some mathematical standard expression that can be calculated through Mathematica or Matlab?

Related bonus question : Can the same method be applied to an asymmetric triangular lattice having dispersion $\epsilon(k_x,k_y)=-2t\,(\cos k_x + \cos k_y)-2t'\,\cos(k_x+k_y)$ ?

• In general, having a larger $0^+$ in the Green's function would reduce the noise. – leongz Mar 4 '17 at 23:27
• Just take more momentum points in your summation, the result can be improved. – Everett You Mar 10 '17 at 9:23

Sorry for the late response but hopefully this can be useful for someone else!

You can reduce the noise using an elliptic integral.

$$D(\varepsilon)=\frac{1}{4 \pi ^3 t}\int_{-\pi }^{\pi }d\phi K\left( \sqrt{1-\left(\frac{\varepsilon +2 t \cos\phi }{4 t}\right)^2}\right)$$

Where K is the Complete Elliptic Integral of the First Kind: http://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html.

It's not trivial to get here. And even from this expression the integral needs to be done numerically with care (it has singularities for many values of $$\varepsilon$$) but it should give better results. Running for six seconds in Mathematica gives me (with $$t=\frac{1}{2\sqrt{3}}$$):

I had the same trouble. I used the formula $$\rho(\epsilon)=\sum_{\vec{k}}\delta(\epsilon - \epsilon(\vec{k}))$$ to numerically calculate the density of states. I did a summation over $$100$$ $$k-$$values for each component, and used a Gaussian distribution with $$\sigma =0.1$$ for the delta function to get the following diagram. Using larger $$\sigma$$ ends up smoothening out the singularity in the derivative at around $$\epsilon=\pm 2$$.

The code was written in C++, and ran for about $$60$$ seconds. $$y$$ axis is $$\rho(\epsilon)$$, $$x$$ axis is $$\epsilon$$ and $$t=1$$.

PS: I did the summation over half the Brillouin zone, which was what I needed for my application.