Quantum versus classical computation of the density of states

Well there is a reason in this case of non-interacting particles- it is the so called "Thermodynamic limit". But I can answer this question without invoking the thermodynamic limit.

One very simple way to see this is using $h$. We know $h ll 1$ so $P = h^{-N}$ for $N gg 1$ will give you $P gg 1$. And in some large limit of the number of particles, you can effectively set $ P rightarrow infty$ i.e. $h rightarrow 0$, which happens to give you the classical limit (This is because setting $h = 0$ gives you a classical theory).

A path integral approach might give a connection. Let's start with the quantum mechanical problem then show how one can take the classical limit. We have $N$, possibly interacting, particles. Let the position of the $i$-th particle be $pmb{x}_i$, we can write the partition function as, $$Z = int prod_i dpmb{x}_i intprod_i Dpmb{x}^prime_i(t) expleft[-frac{1}{hbar}int_0^{betahbar} sum_i left[frac{m}{2} dot{pmb{x}}^{prime 2}_i + V_{text{ext}}(pmb{x}^prime_i) right] + sum_{i>j}V(|pmb{x}^prime_i - pmb{x}^prime_j|) dtau right].$$ Here $dpmb{x}_i = dx_i dy_idz_i$, and $Dpmb{x}_i(t)$ is the path integral measure. In this integral, all paths start and end at the same position, that is $pmb{x}_i(0) = pmb{x}_i(betahbar)$, for all paths.

In the classical limit $betahbar$ is small. In this limit, the particle don't get enough time to move very far away form where they started in the path integral. The reason for this is the kinetic energy term in the Lagrangian, this term make paths with very high velocities not contribute much to the integral, and for small $betahbar$, paths with not very high velocities don't move very far from their staring point. With this approximation, we can treat the potientials as being constant, and take them out of the integral. $$Z = int prod_i dpmb{x}_i expleft[-beta left(sum_i V_{text{ext}}(pmb{x}_i) + sum_{i>j}V(|pmb{x}_i - pmb{x}_j|) right) right] intprod_i Dpmb{x}^prime_i(t) expleft[-frac{1}{hbar}int_0^{betahbar} sum_i frac{m}{2} dot{pmb{x}}^{prime 2}_i dtau right].$$

What is left of the path integral is just a a bunch of free particles. I will spare you the details and use, $$int Dpmb{x}(t) expleft[-frac{1}{hbar}int_0^{betahbar} sum_i frac{m}{2} dot{pmb{x}}^2 dtau right] = left(frac{m}{2pi hbar^2 beta}right)^{3/2},$$ for paths that start and end at the same position. The partition function then become, $$Z = left(frac{m}{2pihbar^2beta}right)^{3N/2} int prod_i dpmb{x}_i expleft[-beta left(sum_i V_{text{ext}}(pmb{x}_i) + sum_{i>j}V(|pmb{x}_i - pmb{x}_j|) right) right].$$

Now I can use the following, $$left(frac{m}{2pihbar^2beta}right)^{3N/2} = frac{1}{h^{3N}}intprod_i dpmb{p}_i expleft[-betafrac{pmb{p}_i^2}{2m} right] $$ to write the partition function as, $$Z = frac{1}{h^{3N}} int prod_i dpmb{p}_i dpmb{x}_i expleft[-beta left(frac{pmb{p}_i^2}{2m} + sum_i V_{text{ext}}(pmb{x}_i) + sum_{i>j}V(|pmb{x}_i - pmb{x}_j|) right) right]. $$ This is nothing but the classical expression for the partition function. In classical mechanics, the partition function is defined up to a constant, however taking this classical limit we know what that constant is. Whenever there is a sum over all states in classical mechanics, we have $frac{1}{h^{3N}} int prod_i dpmb{p}_i dpmb{x}_i$. You can read more about this from Feynman and Hibbs' book Path integrals in Quantum Mechanics, they have a section on statistical mechanics.

For particles in a 1-D potential $V(x)$, we can connect the volume of phase space to the number of quantum states via the WKB approximation. Under the usual WKB assumptions, it can be shown (see, e.g., Liboff or Griffiths) that to have a well-defined wavefunction we must have $$ int_{x_1}^{x_2} p(x,E) , dx = left( n + frac{1}{2} right) frac{h}{2}, $$ where $n$ is an integer, $p(x,E) = sqrt{ 2m(E - V(x))}$, and $x_1$ and $x_2$ are the classical turning points for energy $E$ (i.e., $V(x_1) = V(x_2) = E$.) The curves $pm p(x,E)$ are of course the curves that the particle would take in phase space classically; and so the area enclosed by the (closed) classical trajectory in phase space must be $(n+ frac{1}{2}) h$ for some integer $n$.

Consider now all the states between $n$ and $n + Delta n$. These states span an energy range between $E$ and $E + Delta E$. The area between these curves is the allowed volume of phase space with energies between $E$ and $E + Delta E$; and this is obviously the difference between the corresponding areas enclosed by the energy-$E$ curve and the energy-$E + Delta E$ curve $$ int_{[E, E + Delta E]} dp , dx = left( n + Delta n + frac{1}{2} right) h - left( n + frac{1}{2} right) h = (Delta n) h. $$ This is the shaded area in the diagram above. But by the WKB quantization conditions, there are simply $Delta n = Omega(E; Delta E)$ states in this energy range. Thus, $$ Omega(E; Delta E) = frac{1}{h} int_{[E, E + Delta E]} dp , dx $$ as expected.

It might be possible to generalize this to systems in higher dimensions, but I am not familiar enough with the higher-dimensional versions of these quantization rules to know for sure.