2
$\begingroup$

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!

$\endgroup$
  • $\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 '18 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 '18 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$ – Adam Robert Denchfield Apr 18 '18 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$ – Adam Robert Denchfield Apr 18 '18 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 '18 at 16:28
0
$\begingroup$

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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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