Unitary Integration: A Numerical Technique Preserving the Structure of the Quantum Liouville Equation

The quantum Liouville equation for an $n$-level atomic system driven by external fields has a nontrivial kinematic structure; the quantities $\mathrm{tr}{\rho }^j$, $j=1,2,\dots ,n$ remain constant in time, independent of the Hamiltonian. These invariants are physically significant; the qualitative character of the solution depends on their existence. A generic numerical method will not, in general, preserve these invariants. We present a numerical technique which evolves the density matrix via unitary transformations thus exactly preserving these invariants to all orders in the time step.