  We derive equations for the time evolution of the reduced density matrix of a
collection of heavy quarks and antiquarks immersed in a quark gluon plasma.
These equations, in their original form, rely on two approximations: the weak
coupling between the heavy quarks and the plasma, the fast response of the
plasma to the perturbation caused by the heavy quarks. An additional
semi-classical approximation is performed. This allows us to recover results
previously obtained for the abelian plasma using the influence functional
formalism. In the case of QCD, specific features of the color dynamics make the
implementation of the semi-classical approximation more involved. We explore
two approximate strategies to solve numerically the resulting equations in the
case of a quark-antiquark pair. One involves Langevin equations with additional
random color forces, the other treats the transition between the singlet and
octet color configurations as collisions in a Boltzmann equation which can be
solved with Monte Carlo techniques.
