  We study exclusive quarkonium production in the dipole picture at
next-to-leading order (NLO) accuracy, using the non-relativistic expansion for
the quarkonium wavefunction. This process offers one of the best ways to obtain
information about gluon distributions at small $x$, in ultraperipheral heavy
ion collisions and in deep inelastic scattering. The quarkonium light cone wave
functions needed in the dipole picture have typically been available only at
tree level, either in phenomenological models or in the nonrelativistic limit.
In this paper, we discuss the compatibility of the dipole approach and the
non-relativistic expansion and compute NLO relativistic corrections to the
quarkonium light-cone wave function in light-cone gauge. Using these
corrections we recover results for the NLO decay width of quarkonium to
$e^{+}e^{-}$ and we check that the non-relativistic expansion is consistent
with ERBL evolution and with B-JIMWLK evolution of the target. The results
presented here will allow computing the exclusive quarkonium production rate at
NLO once the one loop photon wave function with massive quarks, currently under
investigation, is known.
