Consider the following passage, via this image:
5.3.1 Density of states
Almost all of the spin-polarized fermionic atoms that have been cooled to ultralow temperatures have been trapped by magnetic fields or focused laser beams. The confining potentials are generally 3D harmonic traps. So let's consider this case in more detail. You might be interested to note that Fermi's original paper on fermionic particles considered this case, not the 3D box case above. As we saw previously, ignoring the zero-point energy in each dimension the eigenvalues (accessible energy states) are given by $\epsilon(n_x, n_y, n_z)=n_x\hbar\omega_x + n_y\hbar\omega_y + n_z\hbar\omega_z$. In order to evaluate the various integrals, we first need to obtain the density of states per unit energy. A rough way to do this is to simply set $k_i=n_i$, so that $$\epsilon^2 = k_x^2(\hbar\omega_x)^2 + k_y^2(\hbar\omega_y)^2 + k_z(\hbar\omega_z)^2 \equiv k^2(\hbar \overline\omega)^2,$$ where $\overline \omega = (\omega_x\omega_y\omega_z)^{1/3}$ is the mean frequency, and $dk_i/\epsilon_i=1/\hbar\overline \omega$. Because $k_i=n_i$ now rather than $k_i=\pi n_i/L$, th 3D density of states is given by $$g(\epsilon) = \frac{k^2}{2} \frac{dk}{d\epsilon} = \frac{\epsilon^2}{2(\hbar\overline\omega)^3}.$$
for the first displayed equation,
shouldn't be $\epsilon^2 =\epsilon_{n_x}^2 +\epsilon_{n_y}^2 + \epsilon_{n_z}^2 + 2\epsilon_{n_x}\epsilon_{n_y} + 2 \epsilon_{n_x}\epsilon_{n_z} + 2\epsilon_{n_y}\epsilon_{n_z}$..?
if I assume $\omega_i=\omega$ for $i=x,y,z$
by
$\epsilon_{n_x}=\hbar \omega n_x $
$\epsilon_{n_y}=\hbar \omega n_y $
$\epsilon_{n_z}=\hbar \omega n_z $
$\epsilon_{n_x,n_y,n_z}=\hbar \omega(n_x +n_y +n_z)$
let $\vec{k}=(k_x,k_y,k_z)$ where $k_i=n_i$
$$\epsilon_{n_x}^2 +\epsilon_{n_y}^2 + \epsilon_{n_z}^2 = \hbar^2 \omega^2 (k_x^2 + k_y^2 +k_z^2 ) = \hbar^2 \omega^2 k^2 \not=\epsilon^2~?$$
And for second displayed equation, why it's not $$\frac{\pi k^2}{2} = \frac{1}{8}4\pi k^2~?$$