# Numerical Sensitivity in Density of States of Tight-binding model

I'm working with the tight-binding model, and I'm trying to learn the basics of how to compute the Density of States (DOS) $N(E)$ numerically.

The DOS is given by

$$N(E) = \frac{1}{N}\sum_k \delta(E-\epsilon_k)\, ,$$

where $\epsilon_k$ is the dispersion relation - in 1D, it has the form $-2t\cos(k)$, where $t$ is a parameter I choose to be 1. We can also write this as an integral,

$$N(E) = \frac{1}{2\pi} \int_k \delta(E - \epsilon(k)) dk$$

We have a couple ways of estimating the DOS. I am taking two approaches, both of which seem highly sensitive to numerical parameters with convergence.

1. Use the Lorentzian approximation to the Delta "function", $\delta(x) \equiv \lim_{\epsilon \rightarrow 0} \frac{\epsilon}{\epsilon^2 + x^2}$. The code is

# Approximate the delta function
def delta(x):
return (1/pi)*(eps/(x**2 + eps**2))


Then simply compute the above sum over my range of $k$s, as shown below for the 1D case.

# Use summation form of density of states for numeric calculation
def N(E):
D = sum([delta(E - disp_e(k1)) for k1 in ks])
# Minimum D for every E should be pi/4 for the 1D case. Unfortunately,   it's going to 0 mostly. Why? How?
dos = (1/Num_sites) * D
return dos

2. Compute the DOS of states as the imaginary part of the Green's function of the system at its poles, $N(E) =-\frac{1}{N\pi} \Im{\sum_k \frac{1}{E - \epsilon_k + i\epsilon}}$, where $\epsilon$ again is supposed to be a small parameter. Here is the code for the functions in 1D, though I've done them for 2D and 3D as well:

# Try to get DOS by summing over k-points of the lattice Green's function
def Greens(E, e_k):
return (1/Num_sites)*(1/(E - e_k + 0.05*1.0j))

def N2(E):
D = (-1/(N*pi))*(sum([Greens(E, disp_e(k1)) for k1 in ks])).imag
return D


Finally, onto my question. I spent quite a while (10 hours) trying to figure out why my code wasn't converging to the 1D analytic solution, despite adding more k-points or having a finer resolution over the energies.

It turned out that despite the definitions, the parameter $\epsilon$ I was choosing was WAY too small. I was choosing numbers on the order of $10^{-4}$, but eventually I figured out that $\epsilon$ only showed correct numerical results when it was around $\epsilon \in [0.1, 0.5]$ for the Green's function method, and between 0.01-0.05 for the delta function method.

Here is the correct solution in 1D, both analytic and numerical, that I finally got. This is for epsilon = 0.01 for the delta function method. Here is what I get when I vary the epsilon parameter to be smaller. The oscillations become divergences at even smaller epsilons and I don't get ANY numerical solution. My question is - why? The definitions imply a small $\epsilon$, and furthermore, I don't see any reason why two separate numerical methods would both not converge for the same reason. Is there stability/convergence analysis for these approximations anywhere?

The reason I want to know is because I want to study systems without analytic solutions, and if the quality of my solution varies erratically with a supposedly arbitrary parameter, I can't be sure of my numerical solution!

• I think that your question might not have caught a lot of attention because you are posing it in a really specific manner. – nicoguaro Apr 17 '18 at 16:16
• Do you have the expression for the analytic result? Also, how are you computing the integral numerically? – nicoguaro Apr 17 '18 at 16:16
• Ah, sorry. I suppose I wasn't sure how to generalize it much. After talking with a professor, I think it comes down to how to integrate nearly singular functions efficiently. – Adam Robert Denchfield Apr 18 '18 at 16:07
• For 1D, the analytic result is basically 1/sin(arccos(x)), where x = E/2t, one of my parameters. Numerically, I am doing basically what I discussed above. Adding up delta functions or discretely adding up the Green's functions. – Adam Robert Denchfield Apr 18 '18 at 16:09
• It would be good if you add those details to the question. I guess that you need to pick the integration method carefully. – nicoguaro Apr 18 '18 at 16:28

Exact Density of states

DosE[x_] := 1/(2 \[Pi]) *1/Sin[ArcCos[-x/2]]


Define delta function

Delta[e0_, \[Sigma]_] :=
1/(Sqrt[2 \[Pi] \[Sigma]^2]) *Exp[-0.5*e0^2/\[Sigma]^2]


let evals be the known eigenvalues of the Hamiltonian matrix, then

DOS[e0_, \[Sigma]0_] := Module[{\[Sigma] = \[Sigma]0, sum},
sum = 0.;
Clear[k];
For[k = 1, k <= Length[evals], k++,
If[Abs[e0 - evals[[k]]] < 6*\[Sigma],
sum = sum + Delta[(e0 - evals[[k]]), \[Sigma]]
];
];
Return[sum/L] (*L=system size*)
]


Plot

Plot[{DosE[x], DOS[x, 0.01]}, {x, -1.99, 1.99}, PlotRange -> All,PlotStyle -> {Black, Red}]


The numerical result (using Mathematica) 100% agree with analytical for L=1024 with 0.01 variance of the Gaussian function (delta function).