Abstract
The Liouville equation governing the evolution of the density matrix for an atomic/molecular system is expressed in terms of a commutator between the density matrix and the Hamiltonian, along with terms that account for decay and redistribution. To find solutions of this equation, it is convenient first to reformulate the Liouville equation by defining a vector corresponding to the elements of the density operator, and determining the corresponding time-evolution matrix. For a system of N energy levels, the size of the evolution matrix is N2 × N2. When N is very large, evaluating the elements of these matrices becomes very cumbersome. We describe a novel algorithm that can produce the evolution matrix in an automated fashion for an arbitrary value of N. As a non-trivial example, we apply this algorithm to a 15-level atomic system used for producing optically controlled polarization rotation. We also point out how such a code can be extended for use in an atomic system with arbitrary number of energy levels.