Quantum corrections in the phase space formulation I'm trying to reconcile the following two statements:

*

*Quantum Mechanics gives physical predictions which are different than the predictions that are obtained in the $\hbar \rightarrow 0$ limit, that is, the predictions of quantum mechanics generally depend on $\hbar$.

*In the phase space formulation of QM, the Moyal product and the quantization map are related by the following formulae:
$$ Q(a) Q(b) = Q(a * b),$$
$$ \text{tr} Q(a) = \frac{1}{2 \pi \hbar} \int dx \int dp a(x, p). $$
Here $Q(a)$ is a quantization of the phase space function $a(x, p)$, and $a * b$ is the Moyal product given by
$$ a * b = a e^{\frac{i \hbar}{2} \omega^{ij} \overleftarrow{\partial_i} \overrightarrow{\partial_j}} b, $$
the arrow notation means that $\overleftarrow{\partial}$ acts on $a$ only and $\overrightarrow{\partial}$ acts on $b$ only.
Now consider a fairly generic problem from classical mechanics. Let's say I have an observable $E(x, p)$, and a probability distribution $\rho(x, p)$ (it can be for example the energy of the system and the Gibbs distribution for some value of $\beta$). The question is: what is the expectation value of $\left< E \right>$? According to classical physics,
$$ \left< E \right> = \int dx \int dp E(x, p) \rho(x, p). $$
Now if we try to compute the quantum corrections to this formula, we end up with
$$ \left< E \right>_{\hbar} = 2 \pi \hbar \text{tr} (Q(E) Q(\rho)) = \int dx \int dp E(x, p) * \rho(x, p). $$
It seems as if the corrections are encoded in the Moyal product. However, that can't be true: it is easy to show that for any $a, b$,
$$ \int dx \int dp \; a * b = \int dx \int dp \; a b. $$
(Demonstrate this by integrating by parts and observing that each non-zero-order term contains a product of antisymmetric $\omega^{i j}$ and symmetric $\partial_i \partial_j b$ and therefore vanishes).
So there are actually no quantum corrections in $\left< E \right>_{\hbar}$ defined above.
On the other hand, it is well known that there are quantum corrections to classical expectation values. Take the expectation value of energy a quantum oscillator for example. Computing it in the Hamiltonian eigenstate basis yields a value $\left<E \right>(\beta, \hbar)$ that depends on $\hbar$ explicitly.
Therefore we have a contradiction, so I must conclude that some of the assumptions I've made before must be incorrect. In particular, I see two possibilities:

*

*The trace formula isn't correct:
$$ \text{tr} ( Q(a) Q(b) ) \neq \frac{1}{2 \pi \hbar} \int dx \int dp \; a * b. $$

*The quantization of the probability distribution isn't $Q(\rho)$ but somehow acquires extra corrections.

Both possibilities seem strange to me.
What's going on here?
 A: (I've done this in "reverse" notation as you in the sense that you have (in my notation) $W_A(x,p)\mapsto \hat A$ so my $W_A(x,p)$ is your $a(x,p)$ and your $Q(a)$ is my $\hat A$.  I suppose it's just "semantics" in a way but I'm also more comfortable in starting with operators and finding the symbols rather than quantizing symbols.)
It is true that, if $\hat A\mapsto W_A(x,p)$, then (up to a constant factor)
$$
\hbox{Tr}(\hat A\hat B)\propto \int dx dp W_A(x,p)\star W_B(x,p) = \int dxdp W_A(x,p)W_B(x,p). \tag{1}
$$
If you already have the symbols, there is no point is using the $\star$-product if you want to do a phase space integration as in (1).
In general, finding the symbols for polynomials is where the action is, although for $xp$ space the $\star$ product has a reasonably manageable form in terms of the exponential as you suggest.  You can then use
$$
 W_{AB}(x,p)=W_A(x,p)\star W_B(x,p)
 $$
to build up symbols for complicated observables from the symbols of the simpler ones.  Of course
$W_A(x,p)\star W_B(x,p)\ne W_B(x,p)\star W_A(x,p)$.
The corrections tend to be more visible in the evolution of an observable or a state since
$$
i\partial_t W_{A}=\{W_H,W_B\}_M:= W_H\star W_A - W_B\star W_H\, .  \tag{2}
$$
For this type of computation, there is no phase space integration and the $\star$-product is unavoidable, even if you already have the symbols.
If $W$ is the Wigner symbol, then the leading term of the Moyal bracket in (2) should basically $i\hbar$ times the Poisson bracket of the symbols $W_H$ and $W_a$, plus the corrections you're looking for,
of size ${\cal O}(\hbar^3)$.  If you use different ordering (say $P$- or $Q$-functions) the corrections are of size ${\cal O}(\hbar^2)$.
A: I figured it out.
Indeed, the quantization of $Q(\rho)$ acquires extra corrections, in a very peculiar way. For example, take the Gibbs distribution for a 1d harmonic oscillator:
$$ \rho(x, p) = e^{- \beta H(x, p)} = e^{- \frac{1}{2} \beta (x^2 + p^2)} $$
(I took $m = \omega = 1$ to simplify the equations).
The crucial observation is that
$$ Q(\rho) = Q(e^{- \beta H}) \neq e^{- \beta Q(H)}. $$
If you take $Q(\rho)$ as your dennsity matrix, like I originally did, you will arrive at the value of $\left< H \right>(\beta)$ (the energy expectation) that doesn't depend on $\hbar$.
The correct density matrix is, however, the operator
$$ e^{- \beta Q(H)} $$
and not
$$ Q(e^{- \beta H}). $$
This is explicitly stated at the start of Wigner's paper.
I still don't understand which physical considerations went into choosing the former expression rather than the latter, however, I now understand where the corrections are coming from.
