Algorithm for solving a pump-probe model for an arbitrary number of energy levels

We describe a generalized algorithm for evaluating the steady-state solution of the density matrix equation of motion, for the pump-probe scheme, when two fields oscillating at different frequencies couple the same set of atomic transitions involving an arbitrary number of energy levels, to an arbitrary order of the harmonics of the pump-probe frequency difference. We developed a numerical approach and a symbolic approach for this algorithm. We have verified that both approaches yield the same result for all cases studied, but require different computation time. The results are further validated by comparing them with the analytical solution of a two-level system to first order. We have also used both models to produce results up to the third order in the pump-probe frequency difference, for two-level systems, and up to first order for three- and four-level systems. In addition, we have used this model to determine accurately, the gain profile for a self-pumped Raman laser, for a system involving 16 Zeeman sublevels in the D1 manifold of Rb87 atoms. We have also used this model to determine the behavior of a single-pumped superluminal laser. In many situations involving the applications of multiple laser fields to atoms with many energy levels, one often makes the approximation that each field couples only one transition, because of the difficulty encountered in accounting for the effect of another field coupling the same transition but with a large detuning. The use of the algorithm presented here would eliminate the need for making such approximations, thus improving the accuracy of numerical calculations for such schemes.