Construct the Density Matrix of a Gaussian State from its First and Second Moments / Wigner Function Borrowing some description for the setup from a question I posted earlier here;
Suppose we have $N$ bosonic modes (or quantum harmonic ocsillators) with the usual commutation relations. Now define the position and momentum operators for the $k^{th}$ mode as, $$\hat{q}_k=\frac{1}{\sqrt{2}}(\hat{a}^\dagger + \hat{a}), \, \hat{p}_k=\frac{-i}{\sqrt{2}}(\hat{a}^\dagger - \hat{a}).$$  Let $\hat{\boldsymbol{X}}=(\hat{q}_1,\hat{p}_1,...,\hat{q}_N,\hat{p}_N)^T,$ then the first statistical moment is the vector of means $\langle \hat{\boldsymbol{X}} \rangle=:\boldsymbol{X}$. The second statistical moment is the covariance matrix $\sigma$ with entries $\sigma_{jk}=\frac{1}{2}\langle \hat{X}_j\hat{X}_k +\hat{X}_k\hat{X}_j\rangle -\langle \hat{X}_j\rangle \langle\hat{X}_k \rangle.$
A Gaussian state is fully characterized by its vector of means and its covariance matrix. And this in fact gives us it its Wigner function as 2N-dim Gaussian distribution.
$$ W_{\boldsymbol{X},\sigma}(\boldsymbol{Y}) =\frac{1}{(2\pi \hbar)^N}\exp\left( \frac{1}{2}(\boldsymbol{Y}-\boldsymbol{X})^T\sigma^{-1}(\boldsymbol{Y}-\boldsymbol{X})\right)\, \text{ Wigner Function for a Gaussian state defined by } \boldsymbol{X},\sigma.$$
A good starting point is the  definition of the wigner function below.
$$ W(\vec{q},\vec{p})=\frac{1}{(2\pi \hbar)^N}\int d^N\vec{x} \,e^{\frac{-i}{\hbar} \vec{p}\cdot\vec{x} }\langle \vec{q} +\frac{1}{2} \vec{x}|\hat{\rho}|\vec{q}-\frac{1}{2}\vec{x}\rangle$$
But I'm not quite sure this will lead to what I'm seeking.
My question is, Given the wigner function $W_{\boldsymbol{X},\sigma}(\boldsymbol{Y})$ mentioned earlier, how can I get the density matrix in the Fock state representation? Explicitly let's assume $N=1$, I want to go from
$$ W_{\boldsymbol{X},\sigma}(\boldsymbol{Y}) \longrightarrow \hat{\rho}=\sum_{m,n=0}^\infty \rho_{mn}|m\rangle\langle n|$$
So how can one obtain $\rho_{mn}$ from the Wigner function.
Any help is appreciated, thanks!
 A: OK, I assume you addressed your question and you are on your way... Just a reminder of language to bridge the distracting communications gap.
For a "Gaussian cigar" (W Schleich's book's nomenclature, not ours!) Wigner function,
$$
W(x,p;a,b;s) =\frac{\pi}{\hbar} \exp \left (    -(x-a)^2/2s^2 -2s^2 (p-b)^2/\hbar^2 \right ),
$$
with $\langle x\rangle =a; \langle p\rangle =b$ and  variance/squeezing of the phase-space ellipse,
σ= diag($s^2,\hbar^2/s^2$), the operator density matrix which produces this Wigner function is given by the Weyl transform of it, eqn (3) of  our booklet,
$$
\hat \rho =\int\!\! dydxdp ~~|x+y/2\rangle ~ W(x,p;a,b;s) e^{ipy/\hbar}~\langle x-y/2| ~,  
$$
in the position representation.
You should plug this into your "starting point definition" to confirm that it, in fact, is but its inverse,
$$
{1\over 2\pi \hbar} \int\!\! dz e^{-iuz/\hbar}\langle q+z/2|\hat \rho|q-z/2\rangle \\
= {1\over 2\pi \hbar} \int\!\! dzdydxdp ~ e^{i(py-uz)/\hbar} W(x,p)\langle q+z/2| x+y/2\rangle \langle x-y/2| q-z/2\rangle\\
={1\over 2\pi \hbar} \int\!\! dzdydxdp ~ e^{i(py-uz)/\hbar}    W(x,p) \delta(q-x)\delta(z-y) \\
={1\over 2\pi  } \int\!\! dw dp ~ e^{iw(p-u) }    W(q,p)  =W(q,u) ~.
$$
You are probably familiar with the standard change of basis from position to energy/number eigenstates,
$$\left\langle x \mid n \right\rangle = {1 \over \sqrt{2^n n!}}~ \pi^{-1/4} \exp(-x^2 / 2)~ H_n(x),$$
as here.
You may gingerly plug in number basis completeness relations on either side of the above density matrix.
Potentially relevant questions here  include  this or that.
