# 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$$?

• How are the $\Psi_n(p)$ stored? – catalogue_number Jun 22 '20 at 13:03
• As list in Wolfram Mathematica: {{$p_1$, $\Psi_n(p_1)$},{$p_2$, $\Psi_n(p_2)$},,..} – Luke Jun 22 '20 at 13:14

## 1 Answer

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$$

• Thank you very much for your response! Is what you are saying equivalent with calculating the sum of my first equation on my momentum grid such that I get a large matrix $\rho_{ab}(t=0)$? I also do not quite understand how you derived the sinc function? – Luke Jun 22 '20 at 16:47
• About your confusions: I tagged the Wigner transform because I thought it might be a way of solving my problem. The $\Psi_n$ are not orthogonal since they arise from different trajectories within the quantum trajectory formalism. My first equation is basically another way of writing: $\rho = \frac{1}{N} \sum_n^N \frac{|\Psi_n\rangle\langle \Psi_n |}{\langle \Psi_n |\Psi_n \rangle }$, where I average over many trajectories. – Luke Jun 22 '20 at 19:57
• The sinc functions arose from $\langle x | \phi_a \rangle = \int dp \langle x | p \rangle \langle p | \phi_a \rangle = \frac{1}{2\pi} \int dp e^{ipx} \text{rect}_{\Delta P}(p-P_a)$ – catalogue_number Jun 23 '20 at 4:51