I would like to visualize $W(X,P)=W(\alpha,\alpha^*)$ from a given density matrix $\hat{\rho}$ that has been obtained before e.g. from the master equation. I am especially interested in density matrices describing a single bosonic mode with Kerr nonlinearity and expressed in Fock states up to some cutoff $N_{\text{cutoff}}\sim 50$. The states are expected to be quite similar to a gaussian state (as a matter of fact, when it is known that the state is truly Gaussian, the Wigner function can readily be obtained from correlation functions, but I explicitly want to see how far this gaussian approximation is from the true state.)

My initial take at calculating the wigner function is using the Characteristic function: $$\chi_W(\lambda,\lambda^*)=Tr\left\{\rho \exp(\lambda a^\dagger-\lambda^*a) \right\}\\ W(\alpha,\alpha^*)=\frac{1}{\pi^2}\int d^2\lambda \exp(-\lambda\alpha^*+\lambda^*\alpha)\,\chi_W(\lambda,\lambda^*) $$ Unfortunately, for let's say a 100x100 grid in the $\alpha$- and $\lambda$-planes the integral sums over 10^8 terms, which takes ages (a day) to evaluate to get a result that is of a too low quality.

I found here a program which uses an alternative approach based on Laguerre polynomials (I haven't found where their algorithm stems from unfortunately). This is already significantly faster (~30 minutes) and gives a meaningful, smooth 'mountain' in phase space corresponding to the state of the system, but aside from that also a giant noisy ring around the origin of fluctuations (I thought for a while this could be physical interference, but it is too irregular and doesn't disappear when feeding a $\hat{\rho}$ that is explicitly gaussian thus classical).

The last source I found is this github page, but it is not entirely transparant to me (partly because I'm unfamiliar with the language). It does tell me that FFT is also an option, though I don't know how to do that for an integral over the whole complex plane.

I prefer to work with MATLAB, but Mathematica is also an option.

  • $\begingroup$ I've always thought Python is reasonably similar to Matlab, but I guess that might be a YMMV type issue? $\endgroup$
    – Kyle Kanos
    Commented Jun 7, 2017 at 20:55
  • $\begingroup$ @KyleKanos I didn't know it was Python (one day, I will learn it). To large extent, I have some idea what is going on, but I don't feel very comfortable distillating an algorithm out of it without derivation. $\endgroup$
    – Wouter
    Commented Jun 9, 2017 at 9:35
  • 1
    $\begingroup$ Well it's changelog you've linked to, the full python file is here $\endgroup$
    – Kyle Kanos
    Commented Jun 9, 2017 at 10:01
  • $\begingroup$ For now, I decided working with the qutip package itself in python, which seems very useful. $\endgroup$
    – Wouter
    Commented Jun 14, 2017 at 12:39

2 Answers 2


I eventually used the Clenshaw-algorithm of the QuTiP toolbox mentioned above and it works adequately and fast for my considered density matrix.


I implement it on mathematica by Laguerre polynomials https://drive.google.com/file/d/1Qg9dW7oS73Y4gUUJCtdONvScDfBY2Rds/view?usp=drivesdk


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.