Consider the following integral in n-dimensions:
\begin{align} g_n(\varepsilon) &= \frac{1}{(2\pi)^n} \int d^n \mathbf k \; \delta(\epsilon +2t\cos (k_xa) + 2t\cos(k_y a)+...) \end{align} The boundaries of this integral is the hypercube with the boundaries $[0,\frac{2\pi}{a})$ along each axis. This integral has a physical meaning but that is not important now. Is it possible to calculate this in $n$-dimensions? Is it even possible to calculate it in two dimension?
I can handle this integral in the 1D case:
$$ g_1(\varepsilon)= \frac{1}{\pi} \int_0^{\pi/a} \delta(\epsilon +2t\cos (ka)) = \frac{1}{2ta\pi} \frac{1}{\sqrt{1-\frac{\varepsilon^2}{4t^2}}} $$
Preliminaries
For $n \in \mathbb{N}$ define the tight-binding energy $E_n \colon \left[0, \frac{2 \pi}{a}\right)^n \to \mathbb{R} , \, E_n (k) = - 2t \sum_{j=1}^n \cos(k_j a) \, .$ Due to the symmetries of the cosine function we may write the density of states as $$ g_n \colon \mathbb{R} \to \overline{\mathbb{R}} \, , \, g_n (\varepsilon) = \int \limits_{\left[0,\frac{2 \pi}{a}\right)^n} \frac{\mathrm{d}^n k}{(2\pi)^n} \, \delta \left(\varepsilon - E_n (k)\right) \stackrel{p = a k}{=} \int \limits_{\left[0,\pi\right)^n} \frac{\mathrm{d}^n p}{2 t (\pi a)^n} \, \delta \left(\frac{\lvert\varepsilon \rvert}{2t} + \sum \limits_{j=1}^n \cos(p_j)\right) \, .$$ Note that $g_n (\varepsilon) = 0$ holds if $\lvert \varepsilon \rvert > 2 t n$, so we will only consider the case $\lvert \varepsilon \rvert \leq 2 t n$ from now on.
The density of states in one dimension
For $n = 1$ a short calculation yields $g_1 (\varepsilon) = \frac{1}{2 \pi a t \sqrt{1 - \left(\frac{\varepsilon}{2t}\right)^2}} $ as mentioned in the question.
The density of states in two dimensions
For $n = 2$ the integral is slightly more complicated but still manageable. We have \begin{align} g_2 (\varepsilon) &= \frac{1}{2 \pi^2 a^2 t} \int \limits_0^\pi \mathrm{d} p_1 \int \limits_0^\pi \mathrm{d} p_2 \, \delta \left(\frac{\lvert \varepsilon \rvert}{2t} + \cos(p_1) + \cos(p_2)\right) \\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\stackrel{u_{1/2} = \pm \cos(p_{1/2})}{=} \frac{1}{2 \pi^2 a^2 t} \int \limits_{-1}^1 \mathrm{d} u_1 \int \limits_{-1}^1 \mathrm{d} u_2 \, \frac{\delta \left(\frac{\lvert \varepsilon \rvert}{2t} + u_1 - u_2\right)}{\sqrt{1-u_1^2} \sqrt{1-u_2^2}} \\ &= \frac{1}{2 \pi^2 a^2 t} \int \limits_{-1}^{1 - \frac{\lvert \varepsilon \rvert}{2t}} \frac{\mathrm{d} u_1}{\sqrt{1-u_1^2} \sqrt{1-\left(u_1 + \frac{\lvert \varepsilon \rvert}{2t}\right)^2}} \\ &\!\!\!\!\!\!\!\stackrel{u_1 + \frac{\lvert \varepsilon \rvert}{4t} = v}{=} \frac{1}{2 \pi^2 a^2 t} \int \limits_{-\left(1 - \frac{\lvert \varepsilon \rvert}{4t}\right)}^{1 - \frac{\lvert \varepsilon \rvert}{4t}} \frac{\mathrm{d} v}{\sqrt{1-\left(v - \frac{\lvert \varepsilon \rvert}{4t}\right)^2} \sqrt{1-\left(v + \frac{\lvert \varepsilon \rvert}{4t}\right)^2}} \\ &= \frac{1}{\pi^2 a^2 t} \int \limits_0^{1 - \frac{\lvert \varepsilon \rvert}{4t}} \frac{\mathrm{d} v}{\sqrt{\left[\left(1 - \frac{\lvert \varepsilon \rvert}{4t}\right)^2 - v^2\right] \left[\left(1 + \frac{\lvert \varepsilon \rvert}{4t}\right)^2 - v^2\right]}} \\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\stackrel{v = \left(1 - \frac{\lvert \varepsilon \rvert}{4t}\right) \sin(\phi)}{=} \frac{1}{\pi^2 a^2 t \left(1 + \frac{\lvert \varepsilon \rvert}{4t}\right)} \int \limits_0^{\pi/2} \frac{\mathrm{d} \phi}{\sqrt{1 - \left(\frac{1 -\frac{\lvert \varepsilon \rvert}{4t}}{1 + \frac{\lvert \varepsilon \rvert}{4t}}\right)^2 \sin^2 (\phi)}} \\ &= \frac{\operatorname{K} \left(\frac{1 -\frac{\lvert \varepsilon \rvert}{4t}}{1 + \frac{\lvert \varepsilon \rvert}{4t}}\right)}{\pi^2 a^2 t \left(1 + \frac{\lvert \varepsilon \rvert}{4t}\right)} = \frac{\operatorname{K} \left(\sqrt{1 - \left(\frac{\varepsilon}{4 t}\right)^2}\right)}{2 \pi^2 a^2 t} \, . \end{align} Here, $\operatorname{K}$ is the complete elliptic integral of the first kind and the last step follows from an identity discussed in this question.
Results in higher dimensions
In three or more dimensions this direct approach does not seem very helpful, as the resulting integrals will be too difficult. We can, however, easily compute the Fourier transform of the density of states in arbitrary dimensions: $$ \hat{g}_n (s) \equiv \int \limits_{\mathbb{R}} \mathrm{d} \varepsilon \, g_n(\varepsilon) \mathrm{e}^{\mathrm{i} s \varepsilon} = \int \limits_{\left[0,\frac{2 \pi}{a}\right)^n} \frac{\mathrm{d}^n k}{(2\pi)^n} \, \mathrm{e}^{\mathrm{i} s E_n (k)} = \left[\int \limits_0^{\frac{2\pi}{a}} \frac{\mathrm{d} k_1 }{2 \pi} \, \mathrm{e}^{- 2 \mathrm{i} t s \cos (k_1 a)} \right]^n = \left(\frac{\operatorname{J}_0 (2ts)}{a}\right)^n \, .$$ Here, $\operatorname{J}_0$ is a Bessel function of the first kind. The inverse transformations for $n \in \{1,2\}$ require some work (they are given in this question), but the results of the direct approach can be reproduced. The case $n = 3$ is the subject of said question and the result is very complicated, so numerical methods are probably better suited to finding the density of states in dimensions $n \geq 3$.