Free evolution of a density matrix in position space I have a density matrix $\rho$ in momentum representation at time $t=0$:
\begin{equation}
  \langle p' |\rho(0) |p\rangle = \sum_{n=1}^{1000} p_n \Psi_n^*(p',0) \Psi_n(p,0) 
      \end{equation}
resulting from a quantum jump calculation (the different $\Psi_n(p)$ are not orthogonal and have large overlaps)
I would like to calculate the time evolution of the position space distribution: $\langle x | \rho(t) |x \rangle$ with the free particle Hamiltonian $H = p^2/4$ without any source of decoherence. $x$ and $p$ are in units of $p_{\mathrm{zp}}$  and $x_\mathrm{zp}$, which are the zero point motions in some harmonic oscillator basis.
I can do that for each individual $\Psi_n(p)$ by calculating
\begin{equation}
\Psi_n(x,t) = \int dp e^{- \mathrm{i} p x/2} e^{- \mathrm{i} p^2/4} \Psi_n(p,0)
\end{equation}
(ignoring pre factors)
and then calculating
\begin{equation}
\langle x | \rho(t) |x \rangle = \sum_{n=1}^{1000} p_n |\Psi_n(x,t)|^2.
\end{equation}
However, doing this takes an extremely long time with my computer as I have the $\Psi_n(p)$ only in numerical form on a discretized momentum space grid.
Is there a clever/much faster way of getting the numerical time evolution of $\langle x | \rho(t) |x \rangle$?
 A: I'm not sure why you've tagged the Wigner transform here - that's only useful for obtaining a quasiprobability x-p distribution.
I think you want the regular old Fourier transform, $$\Psi(x,t) = \frac{1}{2\pi} \int dp e^{ipx} \Psi(p,t)$$
As you probably know, the evolution of the (abstract) density operator is given by the Heisenberg equation,
$$ i \hbar \frac{\partial \rho}{\partial t} = [H, \rho] $$.
To work with this practically, we need to represent it in a basis - For the momentum space representation you have, that basis is essentially rectangular functions corresponding to the grid spacing $\Delta P$, $\phi_j(p) = \text{rect}_{\Delta P}(p-P_j)$ which are obviously orthogonal, but only approximately complete. (the normalisation is also a bit off, but that's just a prefactor) In this sense, you can now reinterpret a state $|\psi\rangle = \sum_j \psi_j |\phi_j\rangle$. Then the initial condition has the matrix representation $$ \rho_{ab}(t=0) = \sum_n w_n \Psi_n (p_a) \Psi^*_n(p_b)$$
for which the N^2-dimensional first-order ODE you must solve is
$$ i \hbar \partial_t \rho_{ab}(t) = [\langle \phi_a | H | \phi_b \rangle, \rho_{ab}(t)]$$
Now, when it comes to actually answering your question - for grid spacing $\Delta P$, grid points P_b, I don't think there's anything faster than this (apart from maybe some trickery with the FFT...
$$ \langle x |\mathbb{1} \rho \mathbb{1}| x' \rangle = \int dp dq \langle x | p \rangle \langle p | \rho | q \rangle \langle q | x' \rangle \\
\approx \sum_a \sum_b \langle x | \phi_a \rangle \langle \phi_a | \rho | \phi_b \rangle \langle \phi_b | x' \rangle \\
= \sum_{ab} \rho_{ab}(t) \frac{e^{iP_a x-iP_b x'}}{\pi \Delta P} \text{sinc}(\Delta P x/2)^2
$$
