Construct the Density Matrix of a Gaussian State from its First and Second Moments / Wigner Function

Question

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!

Answer 1

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.