  We study the production of heavy quarkonium states in high energy
proton-nucleus collisions. Following earlier work of Blaizot, Fujii, Gelis, and
Venugopalan, we systematically include both small $x$ evolution and multiple
scattering effects on heavy quark pair production within the Color Glass
Condensate (CGC) framework. We obtain for the first time expressions in the
Non-Relativistic QCD (NRQCD) factorization formalism for heavy quarkonium
differential cross sections as a function of transverse momentum and rapidity.
We observe that the production of color singlet heavy quark pairs is sensitive
to both "quadrupole" and "dipole" Wilson line correlators, whose energy
evolution is described by the Balitsky-JIMWLK equations. In contrast, the color
octet channel is sensitive to dipole correlators alone. In a quasi-classical
approximation, our results for the color singlet channel reduce to those of
Dominguez et. al. [1]. We compare our results to those obtained combining the
CGC with the color evaporation model and point to qualitative differences in
the two approaches.
