Wigner transform of a thermal density matrix (Schrödinger Hamiltonian) I am interested in computing the Wigner transform in $R^{2n}$ of $$e^{-\beta H}$$
where
$$H = \sum_{k=1}^n P_k^2 + V(x_1, \ldots,x_n)$$
I am assuming that $V$ is a polynomial bounded below  and obviously $$P_k= -i\hbar \frac{\partial}{\partial x^k}$$
I am not necessarily interested in the exact result.  In particular, I would like to know if the Wigner function tends to the corresponding classical exponential for $\hbar \to 0$. Some references on these issues would be welcome.
 A: Your question completely outranges my competence, but Balescu's famous book might be addressing it, at some level. In your shoes, I might choose an oscillator potential and resolve your hamiltonian to decoupled normal modes, which should be tractable with the celebrated (by dint of its non-triviality!) result for just one oscillator, I thought I might review below. It would only be useless if it were universally in mind.
Work in natural units of the problem, $m=1, ~~\omega=1$, i.e. absorb $\sqrt{m\omega}$ into x and 1/p, and $1/\omega$ into H,
rendering it as $H=(p^2+x^2)/2$. The quantum hamiltonian is then $\hat H=(\hat p^2+\hat x^2)/2$, whose Wigner transform is $H=(p^2+x^2)/2$.
The Groenewold-Moyal correspondence of QM operators to phase space variables x and p dictates that the Weyl transform of $\exp(-\beta \hat H)$ is $\exp_\star (-\beta H)$, the star exponential representing formally the power expansion of the exponential with powers supplanted by $\star$-powers.
$$                
\exp_\star \left ( -\beta   H \right )
= \exp_\star \left ( -\beta H \right ) \star 1 \nonumber \\ 
=\exp_\star \left ( -\beta H \right )\star 2\pi\hbar \sum_n f_{n} 
=2\pi\hbar \sum_n e^{-\beta E_n  }f_{n} ~,   \label{spectral}
$$
which is thus a generating function for the stargenfunctions $f_n$.
Of course, for $\beta =0$, the obvious  identity resolution  is recovered.
This Wigner transform, then, is directly seen (Bartlett-Moyal) to sum to
\begin{equation}                
\exp_\star \left (-\beta  H  \right )= 
\left ( \cosh\left(\frac{ \beta \hbar}{2}\right)\right ) ^{-1}
\exp\left ( \frac{-2}{\hbar} \tanh\left(\frac{\beta\hbar}{2}\right) H\right ) ,  
\end{equation}
which is to say just a Gaussian
in phase phase,  with the proper classical $\hbar\to 0$ limit, the starting point.
Details and clarifications could be found in our short booklet, section 0.14 (time evolution).
