  We demonstrate the applicability of a recently proposed multiscale
thermalization algorithm to two-color quantum chromodynamics (QCD) with two
mass-degenerate fermion flavors. The algorithm involves refining an ensemble of
gauge configurations that had been generated using a renormalization group (RG)
matched coarse action, thereby producing a fine ensemble that is close to the
thermalized distribution of a target fine action; the refined ensemble is
subsequently rethermalized using conventional algorithms. Although the
generalization of this algorithm from pure Yang-Mills theory to QCD with
dynamical fermions is straightforward, we find that in the latter case, the
method is susceptible to numerical instabilities during the initial stages of
rethermalization when using the hybrid Monte Carlo algorithm. We find that
these instabilities arise from large fermion forces in the evolution, which are
attributed to an accumulation of spurious near-zero modes of the Dirac
operator. We propose a simple strategy for curing this problem, and demonstrate
that rapid thermalization--as probed by a variety of gluonic and fermionic
operators--is possible with the use of this solution. In addition, we study the
sensitivity of rethermalization rates to the RG matching of the coarse and fine
actions, and identify effective matching conditions based on a variety of
measured scales.
