Do we lose information about the state when we obtain the Wigner function by solving the eigenvalue equation? It can be shown that
$$H(q,p)\star W_{\psi}(q,p)=EW_{\psi}(q,p)$$
where $H(q,p)$ is the classicaly Hamiltonian function, $\star$ is the Moyal/Groenewold star product and $W_{\psi}(q,p)$ is the Wigner function of state $|\psi\rangle$ which is an eigenstate of the Hamiltonian with eigenvalue $E$
$$\hat{H}|\psi\rangle=E|\psi\rangle.$$
This equation can be solved to obtain a PDE whereby the Wigner function is a solutions to this PDE. When applying this to the dimensionless Harmonic oscillator hamiltonian,
$$H=\frac{1}{2}(p^{2}+q^{2})$$
(that is, $m=\omega=1$) it is possible [1] to show that
$$f_{n}(q,p)=e^{-\frac{2H}{\hbar}}L_{n}\left(\frac{4H}{\hbar}\right)$$
are solutions to the PDE. We know from Groenewold's 1949 paper that the Wigner functions  for number states $|n\rangle$ are
$$W_{|n\rangle}(q,p)=\frac{(-1)^{n}}{\pi\hbar}e^{-\frac{2H}{\hbar}}L_{n}\left(\frac{4H}{\hbar}\right).$$  Where $L_{n}$ are Laguerre Polynomials.
In the paper he used the method of inserting the wavefunction into the formula for the Wigner function after solving the Schrodinger equation:
$$W_{|n\rangle}(q,p)=\int_{\mathbb{R}}\frac{dy}{2\pi}e^{-iyp}\langle q+\frac{1}{2}\hbar y|n\rangle\langle n|q-\frac{1}{2}\hbar y\rangle.$$
In principle, I can see how you would obtain the normalisation factor $\frac{1}{\pi\hbar}$ for $f_{n}$ by integrating and demanding that
$$\int_{\mathbb{R}^{2}}W_{|n\rangle}dqdp=1,$$ however I don't see how you could obtain the $(-1)^{n}$ by solving the $\star$genvalue equation, since the functions
$$f_{n}(q,p)= \frac{1}{\pi\hbar} e^{-\frac{2H}{\hbar}}L_{n}\left(\frac{4H}{\hbar}\right)$$
are still solutions to the PDE without this factor. So am I missing something? Or do we generally lose some information about the state when trying to solve for the Wigner function in this way?
 A: In a linear equation in W, any multiple of a solution is also a solution.   This also holds for Schroedinger’s equation, in Hilbert space. Just as in Hilbert space, one needs normalize the solution, in this case, in phase space: $\int dx dp f(x,p)=1$. This normalization is automatic when normalized wavefunctions are input as in Groenewold's  (Wigner's) construction.
In the self-standing solution of the stargenvalue equation, however, one must impose the normalization condition from first principles, since the fundamental conceit of deformation quantization is  that  one has never even heard of a wavefunction or a Schroedinger eqn!
Which is what has been done in equation (65) of the booklet you link.
To see how this works, absorb $\sqrt \hbar$ into the normalization of x as well as p , to see that
$$
f_0= {e^{-(x^2+p^2)}\over \pi} , \qquad f_1=- {e^{-(x^2+p^2)}\over \pi} (1-2(x^2+p^2)), ~~~\leadsto  \\
\int dx dp ~~f_0=1 , \qquad \int dx dp ~~f_1=-(1-2)=1,
$$
etc. for all n. Moments of a 2d Gaussian.  That is, the alternating signs are needed to ensure all WFs are normalized. This is a Sheffer feature of the structure of Laguerre polynomials. Note this is the alternating sign expression of (65), equivalent to your $W_n$ and not your expression. The quickest mnemonic handle on it, as you may check from the figures, is that the hump at the origin is positive for even indices and negative for odd ones.
