I prefer to see it in the following way:
\begin{equation}
Z_{quantum} \equiv \sum_{m} e^{-\beta E_m}
\end{equation}
Where $m$ is a quantum microstate eigenstate of the Hamiltonian. Now, you can split the sum into two parts; a sum over quantum microstates that yields the same energy eigenvalue $E_n$ and a sum over all possible values of $E_n$:
\begin{equation}
Z_{quantum} = \sum_{n = 0}^{\infty} e^{-\beta E_n} \sum_{m \rightarrow E_n} 1 = \sum_{n = 0}^{\infty} g(E_n) e^{-\beta E_n}
\end{equation}
where $n$ labels the different levels of the energy spectrum of the system and $g(E_n)$ is the degeneracy of a given level $n$. When the energy spectrum is evenly spaced, it is quite easy to see what happens: $g(E_n)$ is an increasing function of $n$ and in the thermodynamic limit we can assume that:
\begin{equation}
g(E_n) = e^{S(E_n)/k_B}
\end{equation}
At the end of the day, the sum looks like:
\begin{equation}
Z_{quantum} = \sum_{n=0}^{\infty} e^{-\beta F(E_n)}
\end{equation}
where $F(E_n) = E_n - k_B T S(E_n)$. Since $E_n$ is increasing with $n$ by construction and $S(E_n)$ is an increasing function of $n$ as well, $F(E_n)$ is ensured to have a minimum at some value $n^*$. At $T=0$, the entropy has no weight and the most probable energy is the ground state $n^* = 0$. As $T$ is increased, $n^*$ shifts towards bigger values of $n$. At very high temperature, $n^* \gg 1$ and one can start to think that some approximation may hold.
The idea consists in saying that the partition function is mostly dominated by those states that are in the "neighbourhood" of $n^*$.
Apart from the fact that it is known that at high energies, quantum systems are well described by classical models (coherent states for harmonic oscillator, classical kinetic energy for a particle in a box, Rhydberg states in atoms etc...) one can try something like:
\begin{equation}
Z_{quatum} \sim \sum_{n \sim n^*} e^{-\beta F(E_n)}
\end{equation}Now, since $F(E_n)$ is at an extremum at $n^*$ then it varies very slowly in its neighbourhood which is a a good thing if want to approximate the sum by an integral. To see this, let us consider the term $e^{-\beta F(E_{n+1})} = g(E_{n+1})e^{-\beta E_{n+1}}$. The degeneracy factor $g(E_n)$ is a number associated to $n$ and could also have been called $g_n$. Let us imagine for a minute that $g(E_n)$ is nothing but a continuous function $\Omega(E)$ evaluated at spectrum values $E_n$. We see that imagining such a thing can be helpful because if the energy spacing $\delta E \ll E_n$ then $g(E_{n+1}) = \Omega(E_n+\Delta E) = \Omega(E_n)+\Omega'(E_n)\delta E + \mathcal{O}(\delta E^2)$. Now, $e^{-\beta E_{n+1}}=e^{-\beta E_n}e^{-\beta \delta E} = e^{-\beta E_n}\left(1-\beta \delta E + \mathcal{O}(\delta E^2) \right)$. At the end of the day:
\begin{equation}
e^{-\beta E_{n+1}} = \Omega(E_n)e^{-\beta E_n}\left(1 + (\beta_n-\beta)\delta E \right) + \mathcal{O}(\delta E^2)
\end{equation}where $\beta_n \equiv \partial S/k_B \partial E_n$ is the inverse temperature associated with a system at energy $E_n$ and where I used the fact that $\Omega'(E_n)=\Omega(E_n)\beta_n $. In the thermodynamic limit, as long as the energy levels are in the neighbourhood of $n^*$, $\beta_n$ is very close to $\beta$ and thus, as I qualitatively said earlier, the free energy $F(E_n)$ varies very very slowly. The idea is then to partition the neighbourhood in $m$ lumps $\{\mathcal{L}_{i}\}_{i=1..m}$ of size $\Delta n$ where $F(E_n)$ does not vary. It gives:
\begin{equation}
Z_{quantum} \approx \sum_{i =1 }^{m} \sum_{n \in \mathcal{L}_i} \Omega(E_n) e^{-\beta E_n} = \sum_{i =1 }^{m} \Omega(E_i) e^{-\beta E_i} \Delta n
\end{equation}Let us define then $\rho(E_i)\Delta E \equiv \Omega(E_i)\Delta n$ it finally yields:
\begin{equation}
Z_{quantum} \approx \sum_{i =1 }^{m} \rho(E_i) e^{-\beta E_i} \Delta E \approx \int_0^{+\infty} dE\:\rho(E)e^{-\beta E}
\end{equation}
As I said in the comments, the continuous function $\Omega(E)$ and therefore $\rho(E)$ can be determined using a classical integral over phase space.
With the above derivation, we see that the criterion to be in the classical limit is $\delta E \ll E_n$ (for distinguishable particles). If the system is separable into independent hamiltonians then, this criterion becomes essentially $\delta E \ll k_B T$.
This is nice but unfortunately, it means that one has to compute the full quantum energy spectrum of one particle in the most compliant case (the separable case).
Recently, I have tried to come up with a heuristic way (in the separable case) of getting validity criteria of the classical regime from classical quantities but knowing that the system is genuinely quantum.
The idea is as follows.
The quantum partition function of a particle in a potential $U(\mathbf{x})$ is:
\begin{equation}
z_{quantum} = {\rm Tr}\left(e^{-\beta (\hat{K}+U(\hat{\mathbf{x}}))} \right)
\end{equation}where $\hat{K}$ is the kinetic energy operator. The quantum character of the partition function is embedded in the non commutativity of position and momentum operator that implies that $[\hat{K},U(\hat{\mathbf{x}})] \neq 0$ is not zero. To see that, one can use as a first approximation Glauber's formula (that is a truncated version of the exact Zassenhaus formula):
\begin{equation}
e^{-\beta (\hat{K}+U(\hat{\mathbf{x}}))} \approx e^{-\beta \hat{K}}e^{-\beta U(\hat{\mathbf{x}})}e^{\frac{\beta^2}{2}[\hat{K},U(\hat{\mathbf{x}})]}
\end{equation}If $\beta^2 [\hat{K},U(\hat{\mathbf{x}})] \rightarrow 0$ then, quantum correlations between momenta and positions are effectively vanishing and we get the classical partition function. To see that, the idea is that trace can be computed on any basis set and we choose the basis $\{| \hat{\mathbf{x}} \rangle \}$ we then have:
\begin{equation}
Tr\left[e^{-\beta \hat{K}} e^{ -\beta U(\hat{\mathbf{x}})} \right] = \int \:d\mathbf{x}\: \langle \mathbf{x} |e^{-\beta \hat{K}} e^{-\beta U(\hat{\mathbf{x}}} |\mathbf{x} \rangle
\end{equation}then insert the identity $\mathbf{1} = \int \:d\mathbf{p}\:|\mathbf{p} \rangle \langle \mathbf{p}|$ in between the two exponential functions that yields for the partition function $z_{quantum}$:
\begin{eqnarray}
z &\approx& \int \:d\mathbf{x}\: \langle \mathbf{x} |e^{-\beta \hat{K}}| \int \:d \mathbf{p}\: | \mathbf{p} \rangle \langle \mathbf{p} |e^{ -\beta U(\hat{\mathbf{x}})} |\mathbf{x} \rangle \nonumber \\
&\approx& \int\:d\mathbf{x} d\mathbf{p}\: |\langle \mathbf{x}| \mathbf{p}\rangle|^2 e^{-\beta K} e^{ -\beta U(\mathbf{x})}
\end{eqnarray}By using the fact that $\langle \mathbf{x}| \mathbf{p}\rangle$ is a plane wave with some constant amplitude $C$, we finally get the classical result:
\begin{equation}
z \approx C^2 \:\int\:d\mathbf{x} d\mathbf{p}\: e^{-\beta( K+U(\mathbf{x}))}
\end{equation}
At the end of the day, one has to show that "typical values" of $\beta^2 [\hat{K},U(\hat{\mathbf{x}})]$ are small for the classical limit to hold. One way to so is to use the Robertson inequality that says that for any quantum state $|\psi \rangle$, one has:
\begin{equation}
\beta^2 \sigma_{\hat{K}}(\psi)\sigma_{U(\hat{\mathbf{x}})}(\psi) \geqslant \frac{1}{2}|\langle \psi |\beta^2[\hat{K},U(\hat{\mathbf{x}})] | \psi \rangle|
\end{equation}The problem is now to estimate quantum fluctuations of the kinetic and potential operators. First, if we look at eigenstates of the hamiltonian then, because the energy is constant, then it means that $ \sigma_{\hat{K}}(\psi) \approx \sigma_{U(\hat{\mathbf{x}})}(\psi)$. We thus need a way to estimate kinetic energy fluctuations and this is where it becomes hand-waving again for now. One way to it, I think is to invoke the time-energy Heisenberg relation $\sigma_{\hat{K}}(\psi) \Delta t \geqslant \hbar$. I think that $\Delta t$ is a feature of the system one is looking at and in particular should be related to the potential $U$. Let us call $l$ the typical length scale characterizing the potential $U$ and $v$ a typical value of its speed. We use the fact now that a typical eigenstate will be around the energy level $n^*$ which most likely imply that $v \sim \sqrt{k_B T/m}$. We thus estimate $\Delta t \sim v/l \sim l \sqrt{m/k_B T}$. The most quantum case (called saturation of Heisenberg inequalities) is when $\sigma_{\hat{K}}(\psi) = \hbar/\Delta t \sim \hbar \sqrt{k_B T/ml^2}$. At the end of the day the full criterion is now:
\begin{equation}
\beta^2 \sigma_{\hat{K}}(\psi)\sigma_{U(\hat{\mathbf{x}})}(\psi) \sim \beta^2 \frac{\hbar^2 k_B T}{ml^2}\sim \frac{\hbar^2 }{k_B T ml^2}\ll 1
\end{equation}In the case of an ideal gas, the distance $l$ is nothing but the size of the box $L$ and one finds $L >> \lambda_{beta}$ where I used the notation of Peter here. This formula applies well to all the cases I know but is not very rigorous I admit and again any advice to render it more rigorous is very much welcome.
Best Answer
One way to justify that factor is to consider non-interacting particles in a cubic box with side length $L$. For simplicity let's assume periodic boundary conditions.
The energy eigenstates are also momentum eigenstates. Up to a normalization constant that in general could depend on $p$ (although it will turn out not to in this case), which we will call $\mathcal{N}_p$, the single-particle momentum eigenstates in the position basis are given by \begin{equation} \langle \vec{r} | \vec{p} \rangle = \mathcal{N}_{p} e^{i \vec{p} \cdot \vec{x}/\hbar} \end{equation} Because we are dealing with a box, only certain discrete values of the momentum are allowed, in order that the wavefunction is single-valued when any of the three Cartesian coordinates lying along edges of the box is increased by $L$. For instance, in the $x$ direction, we have \begin{equation} p_x = \frac{2\pi \hbar n_x}{L}, \ \ n_x=1, 2, 3, \cdots \end{equation} To fix the normalization constant, we normalize the states using a Kronecker delta normalization (the advantage of using a box is that we don't have to use a delta function normalization, which resolves some singular behavior). If we label the states by $n_x, n_y, n_z$ (in a hopefully obvious way), then we impose \begin{equation} \langle n_x, n_y, n_z | n_x', n_y', n_z' \rangle = \delta_{n_x, n_x'} \delta_{n_y, n_y'} \delta_{n_z, n_z'} \end{equation} Working through the math leads to the result \begin{equation} \mathcal{N}_p = \sqrt{\frac{1}{L^3}} = \sqrt{\frac{1}{V}} \end{equation} where $V=L^3$ is the volume of the cubic box.
Since the particles are non-interacting, an energy eigenbasis for an $N$ particle state is simply given by a product of $N$ one particle states, thus \begin{equation} \langle \vec{r}_1, \vec{r}_2 \cdots \vec{r}_N | \vec{p}_1 \vec{p}_2 \cdots \vec{p}_N \rangle = \left(\frac{1}{V}\right)^{N/2} \exp\left(i \sum_{a=1}^N \frac{\vec{p}_a \cdot \vec{r}_a}{\hbar}\right) \end{equation} The norm squared of this state is then \begin{equation} |\langle \vec{r}_1, \vec{r}_2 \cdots \vec{r}_N | \vec{p}_1 \vec{p}_2 \cdots \vec{p}_N \rangle|^2 = \frac{1}{V^N} \end{equation} which is the result your book quotes (with slightly more explicit notation on the left hand side).
More generally, this normalization (at least up to a dimensionless constant) follows on dimensional grounds. The normalization of the phase space distribution $p(\vec{r}_1, \cdots, \vec{r}_N, \vec{p}_1, \cdots, \vec{p}_N)$ has to have units of $1/V^N$. For non-interacting particles, physical volume the particles occupy is the only relevant volume scale.