Is there a way to obtain the classical partition function from the quantum partition function in the limit $h \rightarrow 0$? One would like to motivate the classical partition function in the following way: in the limit that the spacing between the energies (generally on the order of $h$) becomes small relative to the energies themselves, one might write
$$Z_{quantum}=\sum_i e^{-\beta E_i} \stackrel{?}{\rightarrow} \int e^{-\beta E}dE$$
(which is clearly wrong on dimensional grounds, but illustrates the idea). Yet this is not of the same form as the classical partition function, namely
$$Z_{classical}=\frac{1}{h^{3N}}\int e^{-\beta H({\bf p},{\bf q})}d{\bf p}d{\bf q}. $$
So one might look for a measure $f(\lambda)$ to append to the quantum partition function such that 
$$\lim_{\lambda \rightarrow 0} \sum_i e^{-\beta E_i} f(\lambda) = \frac{1}{h^{3N}}\int e^{-\beta H({\bf p},{\bf q})}d{\bf p}d{\bf q}.$$
Does such a measure exist, or is there an alternate procedure for deriving the classical partition function from the quantum one?
 A: I assume the system to be $N$ particles with potential energy $U(\vec{R})$ and kinetic energy $\frac{1}{2}(\dot{\vec{R}},M\dot{\vec{R}})$ where $R$ is a $3N$-dimensional vector in the configuration space and $M$ is the mass matrix. In particular, I assume that the particle numbers are fixed -- there is no interactions that create new particles etc. To avoid confusions -- I do consider $m$-particle interactions $U_m(r_1,..,r_m)$, it is all neatly encoded in the function $U$.
Such a system is described by the Lagrangian (note: here I write $+U$ because I need the Wick-rotated Lagrangian there below, general point in quantum statistics)
$$
L=\frac{1}{2}(\dot{\vec{R}},M\dot{\vec{R}})+U(\vec{R})
$$
and the corresponding Hamiltonian
$$
H=\frac{1}{2}({\vec{P}},M^{-1}{\vec{P}})+U(\vec{R})
$$
The quantum partition function is, pretty much by definition
$$
Z=Tr\{\exp(-\beta H)\}
$$
We know that classical limit is conveniently described by saddle points of a functional integral, so it is usefull to reformulate the above formula via a functional integral. This indeed can be done and the result is:
$$
Z=\frac{1}{N!}\int d^{3N}\vec{R}_0\sum_{\sigma\in S_N}\int_{\vec{R}_0}^{\sigma\vec{R}_0}\mathcal{D}[\vec{R}]\exp\left\{-\frac{1}{\hbar}\int_0^{t_\beta}Ldt\right\}
$$
First of all, what are the assumtions. The main assumption is that all particles are identical spinless bosons, there are minor modifications for the case of several species. The case of non-relativistic fermions is also quite similar. The functional integral is taken over $\mathcal{D}[R]$ where particles start at positions $R_0$ at time $t=0$ and end at positions $\sigma R_0$ at time $t=t_\beta=\beta\hbar$. The functional integral itself is pretty much the kernel of $\exp(-\beta H)$ with arguments $R_0$,$\sigma R_0$. The remaining sum and ordinary integral over $R_0$ correspond to taking the trace in the space of bose-symmetric states. Sum over $\sigma$ is the sum over permutations of $N$ particles (fermionic case includes additional factor of $-1$ for odd permutations). Now, the integral is taken over the $R_0$ -- initial positions of the particles.
This integral means the following: take a functional integral where particles leave their initial positions and return back in a short time $t_\beta$. Higher the temperature, shorter the time. To inlude the exchange effects, let the particles exchange their positions. Then take the integral over initial positions. Then you will get the exact quantum partition function.
Now lets turn our intuition on. The time $t_\beta$ contains the Plank constant -- it is 'quantum' small. So, particles in our functional integral that travel to classicaly significant distances have to travel at very high speed in order to return back in time. It means that such trajectories will have extremely large action, so they are supressed  (remeber: we integrate $\exp(-S/\hbar)$). So we have to consider only the trajectories that do not get far away from the initial position $R_0$. In fact, you can show that the characteristic distance from $R_0$ is about the thermal de Broglie wavelength $\lambda_\beta$. We impose the following requirement for the classical approximation: $U$ changes insignificantly on such scales. We also assume that statistically the distances between particles are much larger than $\lambda_\beta$. The latter assumption gives exponential supression to the trajectories that exchange particles. So effectively we only have to consider $\sigma=1$ the identity permutation.
The assumption about $U$ allows us to write, as soon as the particles do not travel very far,
$$
U(\vec{R})\simeq U(\vec{R}_0)+(\nabla U(\vec{R}_0),\vec{R}-\vec{R}_0)+...
$$
the classical approximation corresponds to taking the first term. Higher-order terms correspond to quantum corrections. Let us use this expansion in the functional integral. Then it is clear that the first term can be extracted from under the functional integral, and the functional integral is left only with higher-order terms:
$$
Z\simeq\frac{1}{N!}\int d^{3N} \vec{R}_0\exp\left\{-\frac{t_\beta}{\hbar}U(\vec{R}_0)\right\}\int_{\vec{R}_0}^{\vec{R}_0}\mathcal{D}[\vec{R}]...
$$
Now, if we turn to the classical approximation and forget about all higher-order terms, the functional integral does not depend on $U$ and is easily evaluated to be some constant $C_q$. We are then left with 
$$
Z\simeq\frac{C_q}{N!}\int d^{3N} \vec{R}_0\exp\left\{-\frac{t_\beta}{\hbar}U(\vec{R}_0)\right\}
$$
which is the classical configuration integral. $C_q$ is actually the result of integrating momenta out in the classical partition function (in this problem the momentum integral is Gaussian).
Actually, we can do better and instead of neglecting the higher-order terms, we could have built the perturbation theory. It turns out to be the same at the Wigner-Kirkwood expansion. However, one should be carefull with the exchange trajectories -- they effectively decrease/increase the volume, because when integral over $R_0$ passes through a configuration-space point where two particles are very close to each other (distance around $\lambda_\beta$), the exchange trajectory becomes significant and either amplifies the functional integral (bosons) -- adds up with the normal trajectory or cancels it (fermions). This effect is not exponentially supressed, it contributes to the $\hbar^3$ order (if I recall it correctly).
One last point -- if you want to recover the momentum integral in the classical formula, you should ask yourself: for what purpose? To know about certain averages of momenta. In quantum language, this corresponds to traces of the form $Tr\{f(p)\exp(-\beta H)\}$. To calculate it, consider the generating functional, where instead of $\exp(-S/\hbar)$ we integrate $\exp(-S/\hbar+(\vec{J},M\dot{\vec{R}}(0))$. Then it can be shown to be equivalent to the classical case via the same procedure as above. (I mean that $\langle f(p)\rangle=f(\partial_J)\langle\exp(J p)\rangle|_{J=0}$ in either case).
A: The only trick here is getting used to how discrete sums are turned into integrals.
Suppose you let energy be a function of momentum $p$ and position $q$. Then you can rewrite the discrete quantum partition function as
$Z_{quantum}=\sum_{p,q}e^{- \beta E(p,q)},$
where the sum is over each of the $N$ positions and $N$ momenta, and the only challenge is how to find the appropriate constants for the continuum limit.
This is easiest if you take the system to be in a box of length $L$ and volume $V=L^3$. For position, you want to normalize to the size of the box, ie
$$\sum_q \rightarrow \frac1V\int d^3q$$
for each particle.
For $k$ notice that the spacing of wave numbers in a box is $2\pi n/L$ in each direction. This tells you that the correspondence here is
$\sum_{k} \rightarrow \frac{V}{(2\pi)^3}\int dk$ for each particle.
Put it all together and you get
$$Z_{quantum}=\sum_{p,q}e^{- \beta E(p,q)} \rightarrow \left(\frac{V}{(2\pi)^3}\frac1V\right)^N \int \int e^{- \beta E(p,q)} dq^{3N} dk^{3N},$$
which when you substitute $p=\hbar k$ for each of the 3N k's and collect factors gives you the standard expression.
Note that no special classical approximation was taken here. In fact, classical statistical mechanics is, at least in my view, a misnomer, since you need to use all sorts of things like the discretization of phase space, Planck's constant, the occasional $N!$ factor to avoid Gibbs' Paradox, etc., that make no sense without quantum physics. When using this to derive something like the ideal gas law, the only real classical assumption you make is that Fermi or Bose statistics can be neglected. (This claim seems to be quite disputed in the comments, I'll note, so I will give the disclaimer that this hinges on my personal and somewhat arbitrary consideration of what is considered a 'classical' limit and what is not).
edit: a bit more on the first continuum limit...
Let's take a 1-d discrete system with M sites. Then $\sum_q e^{-\beta E}$ is better written as
$\sum_{i=1}^M e^{-\beta E_i}$, which sums the exponential of energy on each site.
Suppose that the distance between the sites is 'a'. Then $L=Ma$. Furthermore,
$\sum_{i=1}^M =\frac1L\sum_i^M a$
You can probably guess what you want to do now- take a->0 while increasing the number of sites such that L is constant. At this point we can rename a as 'dx' and replace our sum with an integral over it for the identification $\sum_q=\frac1L\int dx$
which when extended to three dimensions and N particles gives the above result.
I certainly won't pretend this is rigorous, but at the same time I think that if you think along these lines you should be able to convince yourself that it couldn't be anything otherwise. Scaling arguments like this come up all over the place, both in statistical mechanics and other areas of physics.
edit2: As Peter rightly points out in the comments, one cannot expand a Hamiltonian simultaneously in the basis of x and p, making it unclear how this classical correspondence should be carried out.
The limit that we are taking is clear enough, I think. In real quantum mechanics, due to noncommutivity each state cannot be thought of as occupying a point in phase space, but rather a probability distribution. In our limit we are assuming that these phase space volumes are small enough to be taken as points- this is another restatement of the continuum limit above.
However, one might reasonably ask for a prescription for how to expand the wavefunction in a basis that treats position and momentum equally, to take this limit. This can be done. The tool used is the Wigner function:
$W_n(x,p)=\frac1h \int_{-\infty}^{\infty} \psi_n^*(x+y) \psi_n(x-y)e^{2ipy/\hbar} dy $
The expectation value of an operator in this formalism is 
$\int \hat{A}(x,p) W(x,p) dx dp$
So if we think of the partition function as $Z_{quantum}=tr(e^{- \beta \hat{H}(p,q)})$
with this formalism in mind and take the limit as before, I think this provides a plausible way to think of relationship between the classical and quantum partition function.
A: 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.
