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

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

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

  • $\begingroup$ 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? $\endgroup$ – Luke Jun 22 '20 at 16:47
  • $\begingroup$ 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. $\endgroup$ – Luke Jun 22 '20 at 19:57
  • 1
    $\begingroup$ 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)$ $\endgroup$ – catalogue_number Jun 23 '20 at 4:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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