  We calculate the yields of quarkonia in heavy ion collisions at RHIC and the
LHC as a function of their transverse momentum. Based upon non-relativistic
quantum chromodynamics, our results include both color-singlet and color-octet
contributions and feed-down effects from excited states. In reactions with
ultra-relativistic nuclei, we focus on the consistent implementation of
dynamically calculated nuclear matter effects, such as coherent power
corrections, cold nuclear matter energy loss, and the Cronin effect in the
initial state. In the final state, we consider radiative energy loss for the
color-octet state and collisional dissociation of quarkonia as they traverse
through the QGP. Theoretical results are presented for J/psi and Upsilon and
compared to experimental data where applicable. At RHIC, a good description of
the high-pT J/psi modification observed in central Cu+Cu and Au+Au collisions
can be achieved within the model uncertainties. We find that measurements of
J/psi yields in proton-nucleus reactions are needed to constrain the magnitude
of cold nuclear matter effects. At the LHC, a good description of the
experimental data can be achieved only in mid-central and peripheral Pb+Pb
collisions. The large five-fold suppression of prompt J/psi in the most central
nuclear reactions may indicate for the first time possible thermal effects at
the level of the quarkonium wavefunction at large transverse momenta.
