von Neumann equation implementation using QuTiP.

QuTiP's mesolve() (master equation solver) integrates the Lindblad master equation. Passing an empty collapse-operator list c_ops=[] disables the dissipator, reducing it exactly to the von Neumann equation:

$$\frac{\mathrm{d}\rho}{\mathrm{d}t} = -\frac{i}{\hbar}[H,\,\rho]$$

Instead of a ket, mesolve() operates on a density matrix rho0 created with qt.ket2dm(). The expectation value of an observable $O$ is now $\langle O\rangle = \mathrm{tr}(\rho\, O)$, computed automatically when O is listed in e_ops. For a pure state with no dissipation, mesolve(c_ops=[]) produces identical trajectories to sesolve().

# run:  python3 von_neumann_qutip.py
# desc: evolve the density matrix rho=|+><+| under H=(omega/2)*Z with no dissipation
 
import numpy as np
import qutip as qt
 
# Initial density matrix: rho = |+><+|
psi0 = (qt.basis(2, 0) + qt.basis(2, 1)).unit()
rho0 = qt.ket2dm(psi0)
 
omega = 2 * np.pi
H = 0.5 * omega * qt.sigmaz()
 
tlist = np.linspace(0, 1.0, 300)
 
# Integrate d/dt rho = -i[H, rho]  (c_ops=[] removes the dissipator entirely)
result = qt.mesolve(H, rho0, tlist, c_ops=[], e_ops=[qt.sigmax(), qt.sigmaz()])
 
# Identical to sesolve() output: <X> completes one full oscillation, <Z> stays 0
print(result.expect[0][0],  result.expect[0][-1])   # ~1.0  ~1.0
print(result.expect[1][0],  result.expect[1][-1])   # ~0.0  ~0.0