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.
In some cases we can approximate the sum by an integral appropriate:
$$\sum_j e^{-\beta\epsilon_j} =\int e^{-\beta E}D(\epsilon)d\epsilon=\int e^{-\beta E(\omega)}D(\omega)d\omega$$
In fact that is how we determine the density function. And I think there is no reason we can write $\sum_i Z_{\omega_i}=\int Z_\omega D(\omega)d\omega$ with the same $D(\omega)$ from the egality above.
And also there is no reason to believe that $Z=\sum_i Z_{\omega_i}$, in fact $Z=\prod_i Z_{\omega_i}$ like you have demonstrated. So that $\ln Z=\sum_i \ln Z_{\omega_i}$, therefore:
$$\ln Z=\int (\ln Z_{\omega_i}) D(\omega)d\omega $$
In fact $D(\omega)d\omega$ counts the number of states having frequency between $\omega$ and $\omega+d\omega$, and similar for $D(\epsilon)d\epsilon$. You can imagine that
$$Z=\prod_{\omega_i}Z_{\omega_i}^{D(\omega)d\omega}$$
which mean we have approximate all the frequency in $[\omega,\omega+d\omega]$ by $\omega$.
Best Answer
In your quantum partition function, you implicitly assumed that the eigenvalues $m^2$ are non-degenerate. In your classical partition function, you have an extra term $$V\int{d^3\vec P\over 2(2\pi)^3\sqrt{m^2+P^2}}$$ that corresponds to the density $g(m)$ of eigenvalue $m^2$ (in the limit of $m$ large). The partition function is $$Z=\int g(m)e^{\beta\lambda m^2}dm$$ Taking into account the degeneracy of the eigenvalues, the quantum partition function should give something similar when going from the sum to the integral.