# 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!

• You neglected to reassure the reader you are completely, unhesitatingly in control of the N=1 case... Label your oscillators by k. Do you appreciate the Weyl map? Commented Mar 16, 2022 at 19:35
• Do not even think of going past N=1, until then. Only then, consult this. Commented Mar 16, 2022 at 19:57
• That's my very point: do not conflate general N with your problem. If, indeed, you can do the N=1 case in your sleep, ask about extension to general N; it is straightforward, as you suspect. But it is not good mental hygiene to throw all pots on the fire without losing track of your problem. Show what you know and have tried for N=1. This is a standard chapter/section in any and all decent textbooks and reviews on the subject. You start with the least usable formulation of the Wigner function, and then roll your eyes that it does not tell you anything about your problem, in your mind. Commented Mar 17, 2022 at 13:45
• Okay, looking at the quoted "position representation" of the Weyl map on the Wiki page and your referenced equation 134 in your book, I now understand what you mean. Yes in this fashion I can use this, and expand the identity operator twice in the fock state representation on the right and left of the density matrix , and integrate the coefficients resulting from the overlap of the position basis states and fock states. That was helpful, many thanks! Though I will have some work to do on Mathematica to do this efficiently :) Commented Mar 17, 2022 at 16:06
• @LostInEuclids5thPostulate a quick word of caution in this arena: sometimes Mathematica runs into challenges with phase space integrals because it doesn't recognize integrals of $\exp(i kx)$ as delta functions, so I always recommend a heavy dose of pen+paper when doing these things alongside Mathematica Commented Mar 17, 2022 at 16:36

## 1 Answer

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.

• Thanks for further posting the details, that's indeed exactly what I did following our exchanges in the comments. I'm now actually working on obtaining the coherent representation (in terms of a complex number alpha as the displacement argument) rather than the typical q,p representation. If you have any pointers on that without spoiling that would be good (I only got started an hour ago, I would like to work on it for a day or two first). Commented Mar 18, 2022 at 18:29
• I'd be shocked if Ulf Leonhardt's book, or Wolfgang Schleich's don't go through this in gory detail... Commented Mar 18, 2022 at 18:31