Van Hove singularity for a two dimensional lattice

Hi everyone. Suppose we consider an electron in a two dimensional lattice, whose dispersion relation is given by:
$$
\epsilon(k_x,k_y)=-J(\cos(k_x a)+\cos(k_y a)),
$$ and where the wave vectors belong to the first Brillouin zone ([itex]k_i\in [-\pi/a,\pi/a][/itex]).

In this case it turns out that the Fermi energy is zero and the Fermi surface is a square defined by:
\begin{align}
&k_y=-k_x+\pi/a\;\;\;\;\text{ for the first quarter of the k-space} \\
&k_y=k_x-\pi/a\;\;\;\;\text{for the second quarter} \\
&k_y=-k_x+\pi/a\;\;\;\;\text{for the third quarter} \\
&k_y=k_x+\pi/a\;\;\;\;\text{for the fourth quarter}
\end{align}
I would like to compute the density of states near the Fermi surface. The density of states is given by:
\begin{align}
\nu (\epsilon)=\int\frac{d\vec l}{4\pi^2}\frac{1}{|\nabla \epsilon|}.
\end{align}
In the situation I am describing the we have, along the Fermi surface:
\begin{align}
|\nabla \epsilon|=|\sin(k_x a)|.
\end{align}
Since I should integrate from [itex]-\pi/a[/itex] to [itex]\pi/a[/itex] the integral is clearly divergent. This is the so called Van Hove divergency. How do I deal with it?

I am expecting a logarithmic divergency around [itex]\epsilon\simeq0[/itex], how do I find it?

Actually both, the divergence for which [itex]\nabla \epsilon=0[/itex] should occur at [itex]\epsilon\simeq\epsilon_F=0[/itex]. I will try the Taylor series, but the problem is that [itex]\epsilon\simeq0[/itex] is realized not only when [itex]k_x\simeq0[/itex], so it seems to me like I'm losing informations.

Actually both, the divergence for which [itex]\nabla \epsilon=0[/itex] should occur at [itex]\epsilon\simeq\epsilon_F=0[/itex]. I will try the Taylor series, but the problem is that [itex]\epsilon\simeq0[/itex] is realized not only when [itex]k_x\simeq0[/itex], so it seems to me like I'm losing informations.

Ah, ok. Yeah, that Taylor series (unless you do the full infinite sum, in which case you wouldn't lose information) is mainly useful near the origin. But the other solutions are for [itex]k_x =\pm \pi/a[/itex], right? And [itex]|\nabla \epsilon|[/itex] is an even function, isn't it?

Yes, in principle you just need to integrate over [itex]k_x\in[0,\pi/a][/itex] since all the other integrals are equal. The point is that in this range the integral is divergent for [itex]k_x=0,\pi/a[/itex]. My question was how to express such a divergence as a function of the energy $\epsilon$ instead of the wave number [itex]k_x[/itex].

Yes, in principle you just need to integrate over [itex]k_x\in[0,\pi/a][/itex] since all the other integrals are equal. The point is that in this range the integral is divergent for [itex]k_x=0,\pi/a[/itex]. My question was how to express such a divergence as a function of the energy $\epsilon$ instead of the wave number [itex]k_x[/itex].

I see. Yeah, I'm not sure about how to do that. But didn't we use [itex]dE=\nabla E \cdot \mathbf{k}[/itex] to remove one of the momentum dimensions already? That kind of suggests that it might be simplest to deal with the integral in terms of the wave numbers. Then again, I might be wrong.

In any case, we should be able to split the integral into different regions, close to each singular point. Then one could do Taylor expansions about the different poles, which shows that all these points have logarithmic divergences and will dominate the integral. (Since there are three of them, their contributions can't cancel out completely.) Then, in a very roundabout way, you could argue that the zero energy case corresponds to the sum of these points, which are associated with logarithmic divergences of k, and thus of the sines and by your argument above, this is related to both the energy gradient and the energy itself. It'd certainly be more elegant if it was possible to rewrite it in terms of an energy integral though.

Yes, this is probably the most proper way to deal with it. My idea is that since I am integrating over the constant energy path [itex]\epsilon (k_x,k_y)=\epsilon[/itex] the energy dependence should come out from the final result of the integration. However I can't see how.