As has been noticed in other answers, the brute force method is quite slow. One way to deal with the evaluation of a large sum of delta functions is to expand the density of states in terms of orthogonal polynomials. Since the range of energies is finite, the Chebyshev polynomials $T_n(x)$ is an obvious choice. Quite generally, a function $f:[-1,1]\rightarrow \mathbb{R}$ can be written as a sum
$$
f(x)=\frac{1}{\pi\sqrt{1-x^2}}\left[\mu_0+2\sum_n\mu_nT_n(x)\right],
$$
where $\mu_n$ are the moments computed according to:
$$
\mu_n=\int_{-1}^1f(x)T_n(x)dx.
$$
It is clear that in the case of $f$ given by a sum of delta functions, the integral is trivial. In principle, this already gives a prescription. However, when the polynomial sum is truncated, the resulting density of states suffers from fluctuations—also known as Gibbs oscillations. One way to deal with this problem is to suppress the weight of higher order polynomials with the $g_n$ factors:
$$
f(x)=\frac{1}{\pi\sqrt{1-x^2}}\left[\mu_0g_0+2\sum_n^N\mu_ng_nT_n(x)\right].
$$
In Rev. Mod. Phys. 78, 275 (2006), it is shown what that the optimal form of the factors is:
$$
g_n=\frac{(N+2-n) \cos \left(\frac{\pi n}{N+2}\right)+\cot \left(\frac{\pi }{N+2}\right) \sin \left(\frac{\pi n}{N+2}\right)}{N+2}.
$$
This leads to the following code:
DOScb[evals_,Ncb_,Nint_]:=Module[{eSorted,neng,xvals,emn, emx,μcb,gcb,itab,de},
eSorted=Sort[evals];
emn=eSorted[[1]]-0.5;
emx=eSorted[[-1]]+0.5;
neng=Length[eSorted];
xvals=2(eSorted-emn)/(emx-emn)-1;
de=2./(Nint-1);
μcb=Table[Total[ChebyshevT[n,#]&/@xvals]/neng,{n,Ncb}];
gcb=Table[N[1/(Ncb+2)*((Ncb+2-n)*Cos[π*n/(Ncb+2)]+Sin[π*n/(Ncb+2)]Cot[π/(Ncb+2)])],{n,Ncb}];
itab=Table[{emn+(x+1.)*(emx-emn)/2,2/(emx-emn)*1/(π*Sqrt[1-x^2])*(1+2*Sum[μcb[[n]]*gcb[[n]]*ChebyshevT[n,x],{n,Ncb}])},{x,-1.+de,1.-de,de}];
Interpolation[itab,"ExtrapolationHandler"->{0&,"WarningMessage"->False}]
]
The function takes 3 arguments: a set of energies, the highest order of included Chebyshev polynomial, and the number of points for interpolation of the final result.
By chance, the exact analytic solution is also known:
$$
n(E)=\frac{1}{2 \pi ^2}K(1-(E/4)^2) \theta (4-|E| ),
$$
where $\theta$ is the Heaviside function and $K$ is an elliptic integral.
Below, I compare density of states computed using 250 Chebyshev polynomials with the analytic expression
t = 1.;
energy[kx_, ky_] := -2 t (Cos[kx] + Cos[ky]);
nkp = 150;
dk = 2*Pi/(nkp - 1);
energyValues =
Flatten[Table[
energy[kx, ky], {kx, -Pi, Pi, dk}, {ky, -Pi, Pi, dk}]];
dos2d = DOScb[energyValues, 250, 1001]
dosExact[w_] := HeavisideTheta[4 - Abs[w]]/(2*Pi^2)*EllipticK[1 - (w/4)^2]
Plot[{dosExact[w], dos2d[w]}, {w, -4.5, 4.5},
PlotTheme -> {"Frame", "Monochrome"},
PlotRange -> All,
Exclusions -> None,
PlotStyle -> {Automatic, Red}]
DiracDeltais the implementation of the $\delta$-distribution in Mathematica, not a usual function. Its plots make no sense. $\endgroup$t? How can you plot anything when you miss the definition oft? $\endgroup$