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. Hubbard U = 0 1D model, epsilon = 0.01.

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. Epsilon = 0.001. At even smaller values the solution diverges.

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!

  • $\begingroup$ I think that your question might not have caught a lot of attention because you are posing it in a really specific manner. $\endgroup$
    – nicoguaro
    Apr 17, 2018 at 16:16
  • $\begingroup$ Do you have the expression for the analytic result? Also, how are you computing the integral numerically? $\endgroup$
    – nicoguaro
    Apr 17, 2018 at 16:16
  • $\begingroup$ 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. $\endgroup$
    – Slenderman
    Apr 18, 2018 at 16:07
  • $\begingroup$ 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. $\endgroup$
    – Slenderman
    Apr 18, 2018 at 16:09
  • $\begingroup$ It would be good if you add those details to the question. I guess that you need to pick the integration method carefully. $\endgroup$
    – nicoguaro
    Apr 18, 2018 at 16:28

2 Answers 2


You can find discussion of this method in the context of particle methods for advection PDEs here. Some of the results at the end of this paper will illuminate the problem. Essentially, as $\epsilon\to0$, your solution will converge in the integral sense, but pointwise estimation of this function is not so easy. If $\epsilon$ is too large, then the solution is oversmoothed and if $\epsilon$ is too small, then you see this jagged behavior.

In particular, there is likely an error estimate for your problem mirroring Theorem 2.2 in the linked paper. Your discretization of the density gets more accurate as you take more points, and the closer these points are, the smaller you are able to take $\epsilon$. This is because for smooth approximations of the Dirac delta function, the parameter $m$ in the paper equals $\infty$. This happens with your Cauchy approximation as well as Gaussian approximations to the Dirac delta function. with $m=\infty$, it is then essential to have $h/\epsilon$ securely less than 1, where $h$ is the separation between your Dirac masses.

A common choice in particle methods is something like $\epsilon = \frac{1}{4}\sqrt{h}$, which ensures that this prefactor is small regardless of choice of $h$ sufficiently small. The underlying analysis is obviously different from your case, which may result in slightly different convergence orders, but the intuition remains the same. You should be able to fix this by playing around with how $\epsilon$ should change along with the number of points you are using.


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.;
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[{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).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.