Appendix B — From classical densities to the Wigner function

So far we have:

• a probability density \(\rho(x,p,t)\) for classical ensembles, obeying Liouville’s equation;
• a wavefunction \(\psi(x,t)\) or a density operator \(\hat\rho\) for quantum systems, obeying the Schrödinger/von Neumann equations.

The Wigner function builds a bridge between these two pictures.
It lives in phase space like \(\rho(x,p)\), yet it is derived from \(\hat\rho\) and keeps all quantum information.

B.1 Definition

For a one–dimensional system the Wigner function is

\[ W(x,p,t) \;=\; \frac{1}{2\pi\hbar}\, \int_{-\infty}^{+\infty} dy\; e^{-i p y/\hbar}\; \bigl\langle x+\tfrac{y}{2}\bigl|\hat\rho(t)\bigr|x-\tfrac{y}{2}\bigr\rangle. \tag{B.1}\]

A quick checklist:

• \(x\) and \(p\) are simultaneous variables (even though they do not commute quantum mechanically).
• \(W\) is real, but it can take negative values – an unmistakable quantum signature.
• The normalisation matches that of \(\rho\): \(\displaystyle\int dx\,dp\,W(x,p,t)=1\).

Classical ↔︎ Quantum analogy

• \(\rho(x,p,t)\) is always non-negative.
• \(W(x,p,t)\) can be negative, revealing quantum interference.

B.2 Marginals and expectation values

Despite possible negativity, \(W\) returns the correct probability densities for position and momentum:

\[ \int_{-\infty}^{+\infty} dp\;W(x,p,t) \;=\; \langle x|\hat\rho|x\rangle \;=\; |\psi(x,t)|^{2}, \]

\[ \int_{-\infty}^{+\infty} dx\;W(x,p,t) \;=\; \langle p|\hat\rho|p\rangle . \]

Expectation values of symmetrised operators follow the phase-space average rule

\[ \langle \hat A \rangle \;=\; \int dx\,dp\;A_{\text{W}}(x,p)\;W(x,p,t), \]

where \(A_{\text{W}}\) is the Weyl symbol of \(\hat A\) (the phase-space version of the operator).

B.3 Properties worth remembering

Property	Classical \(\rho\)	Quantum \(W\)
Real	✓	✓
Non-negative	✓	✗ possible negativity
Normalised	✓	✓
Obeys continuity / Liouville	✓	✓, plus \(\hbar\)-corrections
Supports interference fringes	✗	✓

B.4 Example : Gaussian wavepacket

For a minimum-uncertainty Gaussian wavepacket

\[ \psi(x) \;=\; \frac{1}{(2\pi\sigma_x^2)^{1/4}} \exp\!\Bigl[-\frac{(x-x_0)^2}{4\sigma_x^2} + i\,p_0(x-x_0)/\hbar\Bigr], \]

the Wigner function is also Gaussian and everywhere positive:

\[ W(x,p) \;=\; \frac{1}{\pi\hbar}\, \exp\!\Bigl[-\frac{(x-x_0)^2}{2\sigma_x^2} -\frac{2\sigma_x^2}{\hbar^2}(p-p_0)^2\Bigr]. \]

Negativity appears only when the state contains quantum interference, for example in superpositions of spatially separated Gaussians.

B.5 Why the Wigner function matters

• It lets us visualise quantum states in the familiar \((x,p)\) plane.
  
• Many semiclassical techniques expand around \(W\) and truncate the series.
  
• In quantum optics the Wigner function of an electromagnetic mode can be measured with homodyne tomography (see later lectures).

Using \(W(x,p,t)\) we now have a complete trio:

1. Classical ensemble — \(\rho(x,p,t)\) (Liouville).
   
2. Quantum wavefunction/density operator — \(\psi(x,t)\) / \(\hat\rho(t)\) (Schrödinger/von Neumann).
   
3. Quantum phase-space picture — \(W(x,p,t)\) (Moyal evolution).

B.6 Example: the nonlinear oscillator

In Chapter 7 we solved the Liouville equation for a classical nonlinear oscillator with Hamiltonian

\[ H = \frac{p^2}{2m} + \frac{1}{2} k x^2 + g x^4. \tag{B.2}\]

Moreover, in Chapter 10 we solved the Schrödinger equation using QuTiP for aq simple harmonic oscillator. We now combine these two approaches to study the quantum nonlinear oscillator with Hamiltonian Equation B.2, but with the quantum operator \(\hat x\) and \(\hat p\) instead of the classical variables \(x\) and \(p\).

import numpy as np
from qutip import *

N = 120
m = 0.5 # Mass of the particle
k = 2.0 # Spring constant
G = 0.15 # Nonlinear constant
w = np.sqrt(k/m) # Angular frequency

a = destroy(N)

x = np.sqrt(m * w / 2) * (a + a.dag())
p = 1j * np.sqrt(m * w / 2) * (a.dag() - a)

H = p**2 / (2 * m) + k * x**2 / 2 + G * x**4

# Initial state: coherent state
alpha = np.sqrt(1 / (2 * m * w)) * 1 + 1j * np.sqrt(m * w / 2) * 0.1
psi_0 = coherent(N, alpha)

tlist = np.linspace(0, 2, 500)
result = sesolve(H, psi_0, tlist)
result

<Result
  Solver: sesolve
  Solver stats:
    method: 'scipy zvode adams'
    init time: 0.00020766258239746094
    preparation time: 0.00016641616821289062
    run time: 0.9459991455078125
    solver: 'Schrodinger Evolution'
  Time interval: [0.0, 2.0] (500 steps)
  Number of e_ops: 0
  States saved.
>

The Wigner function can be computed from the resulting state using QuTiP’s wigner function. Let’s plot its evolution over time:

import matplotlib.pyplot as plt
from IPython.display import HTML
from matplotlib.animation import FuncAnimation

# Phase space grid
x_list = np.linspace(-5, 5, 120)
px_list = np.linspace(-5, 5, 120)

plt.rcParams.update({"font.size": 8})
fig, ax = plt.subplots(figsize=(4.6, 2.8))

fig.suptitle(r"Wigner function evolution of a nonlinear oscillator")

img = ax.pcolormesh(x_list, px_list, wigner(result.states[0], x_list, px_list),
                    shading="gouraud", rasterized=True,
                    vmin=-0.3, vmax=0.3, cmap="PuOr")
ax.set_xlabel("Position $x$")
ax.set_ylabel("Momentum $p_x$")

def animate(i):
  wig_i = wigner(result.states[i], x_list, px_list)
  img.set_array(wig_i.ravel())
  return img,

ani = FuncAnimation(fig, animate, frames=range(0, len(tlist), 10), interval=200, blit=True)

plt.close(fig)
HTML(ani.to_jshtml())

Once Loop Reflect