Efficient numerical evaluation of Wigner function

Question

Suppose we want to calculate the Wigner function of some state $|\Psi\rangle = \sum_{n=0}^{N_{max}} c_n|n\rangle$ ($|n\rangle$ are the eigenstates of the Harmonic oscillator) numerically. Starting from the definition,

$$ W(x,p) =\frac{1}{2\pi}\int_{-\infty}^\infty e^{-ip\xi}\Psi^*(x-\xi/2)\Psi(x+\xi/2) $$

I have seen that when the states $\Psi$ which in my case are discrete, people use the Fast Fourier Transform (instead of evaluating directly the integral). Could someone describe how to carry $W(x,p)$ to an evaluation of a discrete Fourier transform in two dimensions?

Ps. Incidentally, I don't want to deal with the explicit evaluation of the integral. I know that the Wigner function for the harmonic oscillator eigenfunction is in terms of a gaussian plus a generalized Laguerre polynomial, but this is of not use to me.

EDIT:

Actually the method that I am asking is already implemented in here; however, it is not clear to me why it is done the way it is.

Answer 1

This is a bit more general since it's for density matrices, and perhaps it is not the fastest way for pure states, but here goes:

Suppose you have a finite difference grid for your spatial coordinate ($x$). In this approach you will first need to express the density matrix in spatial coordinates, i.e., $\rho_{xy} = \langle x|\rho|y\rangle$.

(If you are starting out with a density matrix in the energy eigenbasis (i.e. you know Hamiltonian is diagonalized as $H = P D P^{-1}$), then you can easily convert that energy-basis density matrix to the spatial-basis with some matrix multiplications by $P$ and its inverse. This may actually be the slowest step, given how slow matrix multiplications are.)

Since you're starting out with a pure state, you just need to express it as a spatial wavefunction, then do the usual self-outer-product thing to get the spatial density matrix. This is fast.

Now, here's a trick for doing Wigner function from that spatial density matrix. If you look carefully at the Wigner function integral, you'll notice it's simply a Fourier transform along the antidiagonals of that matrix. To do Fourier transform along antidiagonals, here's a good way:

1. Convert the antidiagonals to columns: Perform a 'shearing' on the matrix, shifting the second row by one element, shifting the third row by two elements, etc.. You'll need to make space for the moved rows, so your matrix is going to end up twice as wide as you started with, i.e., $N \times N \rightarrow N \times (2N-1)$.
2. Perform the discrete Fourier transform along the columns.
3. Pixel shift: As usual for DFT you need to recenter the result, in this case, rolling each column around by half.
4. Phase shift: Now you need to perform a complex phase shifting of each element to compensate for the fact that the zero of each column (the diagonal element from the original matrix) was not centered. So multiply each matrix element by an appropriate phase factor, which effectively performs a vertical shift in the un-fft'd matrix. You could also have just rolled the columns first but for the odd antidiagonals, that would require a half-pixel shift (easier to do after DFT!).
5. Merge even/odd columns: there were actually two kinds of antidiagonals, some which contain diagonal elements, some which don't. Both contain meaningful information, however if your grid is dense enough (your space-basis density matrix was smooth) then they will have very similar information. You can average them together or throw out the weird ones, whatever you like.

I used all of the aforementioned tricks in the following animation. I've included a still frame from it below. The source code is in the link. I have no idea if it is standard or the best way, but I like it. :D

** https://commons.wikimedia.org/wiki/File:Hamiltonian_flow_quantum.webm **