I'm interested in simulating the time evolution of a Wigner function for a harmonic oscillator (and possibly some other potentials) and I can't seem to find a good resource for that. My background in numerical methods is rather limited, so I'm looking for something suitable for a beginner (explanation of algorithms used, for instance).


By the Wigner function I mean the quasiprobability distribution used for quantum mechanics in phase space, as defined here http://en.wikipedia.org/wiki/Wigner_quasiprobability_distribution

The time evolution equation I'm interested in solving can be expressed in many different ways, though in the case of the harmonic oscillator it reduces to the Liouville equation:

$$\frac{\partial W(x,p;t)}{\partial t} = -\frac{p}{m}\frac{\partial W(x,p;t)}{\partial x}+m\omega^2x\frac{\partial W(x,p;t)}{\partial p}$$

This is equivalent to the following:

$$\frac{\partial W(x,p;t)}{\partial t} = -\frac{p}{m}\frac{\partial W(x,p;t)}{\partial x}-\frac{1}{2 \pi \hbar^2} \int_{-\infty}^{\infty}dp'F(x,p-p')W(x,p')$$

with the kernel

$$F(x,p-p')=\int_{-\infty}^{\infty} dx' \sin \left( \frac{x'(p-p')}{\hbar} \right) \left[ V \left( x+\frac{x'}{2} \right)- V \left( x-\frac{x'}{2} \right) \right]$$

where $V(x)=\frac{1}{2}m \omega^2x^2$

So that $$F(x,p-p') =m\omega ^2 x\int_{-\infty}^{\infty} dx' x' \sin \left( \frac{x'(p-p')}{\hbar} \right). $$

Although this doesn't converge on its own.

  • 2
    $\begingroup$ It may be helpful to provide the reader with what a Wigner function is, as nomenclature for functions sometimes vary across timezones. Generally, most programs deal with solving some sort of differential equation, so providing that will be helpful as well. $\endgroup$
    – John M
    Nov 1, 2014 at 15:30
  • $\begingroup$ Ok sure, I did that. I guess for the simplest case (harm. oscillator) solving the Liouville is probably the easiest, though I'm interested in learning how to deal with the integral equation in practice. $\endgroup$ Nov 1, 2014 at 16:02
  • $\begingroup$ The lack of convergence definitely seems to be an issue. Where did you find such a kernel? $\endgroup$
    – Kyle Kanos
    Nov 5, 2014 at 15:03
  • $\begingroup$ Unless I made some sort of simple calculational error, this kernel is equivalent to the series representation: (44b) here stat.physik.uni-potsdam.de/~pikovsky/teaching/stud_seminar/… . Writing the sine as an exponential and substituting the definition of the Wigner function leads to that result. $\endgroup$ Nov 5, 2014 at 16:14

3 Answers 3


Your equation in the Liouville form is elementary for numerical integration, it is structurally just a linear advection equation with spatially varying coefficients. The transformed equation with the kernel F is not useful at all for numerical solution, don't bother with it.

All we have here is a 2D advection equation (I use y instead of p):

$$ \partial_{t} W(x,y,t) = V_x(x,y) \partial_{x} W(x,p,t) + V_y(x,y) \partial_{x} W(x,p,t) $$

where $V_x=C_1 y$ and $V_y=C_2 x$.

A simple way to solve it is using explicit time-stepping for the spatially and temporally discretized form

$$ (W_{i,j}^{k+1}-W_{i,j}^{k})/\tau = V_x (W_{i,j}^{k}-W_{i-1,j}^{k})/\delta x + V_x (W_{i,j}^{k}-W_{i,j-1}^{k})/\delta y $$

This elementary time-stepping scheme will allow time integration; as long as you provide boundary conditions upstream you can march along the characteristics. Note that in this method some constraints on the time step will arise from the stability requirements, and the accuracy order is low. Also note that the upwind spatial discretization is used for the advection operators, otherwise it will be unstable for any time step. You can improve this a bit by using explicit Lax-Wendroff or some other classical scheme for hyperbolic PDEs. However, with today's computers, you can as easily use implicit time-stepping and just solve for the whole domain as a linear system of equations. To improve the accuracy I would go to the 4th order in the spatial discretization. Note that $(x=0,y=0)$ is a singular point here so some care should be taken to avoid using it in the numerical stencil.

  • $\begingroup$ Thank you very much for your answer. However, since I have next to no experience in numerical methods, I don't understand a few things. First of all, what is meant by "as long as you provide boundary conditions downstream you can march along the characteristics" and "4th order in the spatial discretization"?? $\endgroup$ Nov 5, 2014 at 22:42
  • 1
    $\begingroup$ It is hard to explain all this here, you just need to read a basic computational book, e.g., Numerical Heat Transfer and Fluid Flow by Patankar, that will give you enough background knowledge for this. $\endgroup$ Nov 6, 2014 at 0:24

Well, for the harmonic oscillator, you have the full closed answer, so you don't really need numerics. It is, as Groenewold discovered in his 1946 thesis (Thomas L. Curtright, David B. Fairlie, & Cosmas K. Zachos, A Concise Treatise on Quantum Mechanics in Phase Space, World Scientific, 2014. The PDF file is available here.), merely rigid rotation!:

$$W(x,p; t) = W(x \cos t − p \sin t, p \cos t + x \sin t; 0).$$

Have set $m=ω=1$ for simplicity.

A state of the art numerical simulation example this year is

The following classic references rely on numerical simulation to draw sound conclusions with useful insights into WF time evolution, and include links to numerical sources:


I don't know if you are still interested in this question, but there is some work in the literature that might be very interesting for you:

In this paper they propose an interpretation of the Wigner function as a particular wave function and in this other paper they propose a numerical method to compute the propagation, note that you might be interested in the older version of the paper (v2) which also contains python code in the appendix!


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.