  NLO evolution of the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner
(JIMWLK) equation with massless quarks was derived a few years ago. We make a
step further to compute the evolution kernels focusing on the effects due to
finite quark masses. To this goal, the light-cone wave function of a fast
moving dilute hadronic projectile is computed up to ${\cal O}(g^3)$ in QCD
coupling constant. Compared with the massless case, a new IR divergence
emerges, which is eventually canceled by a mass dependent counter term. Our
results extend the theoretical tools used in physics of gluon saturation and
aim at improving precision in future phenomenological applications.
