  Numerical simulations of the quarkonium spin splittings are done in the
framework of lattice nonrelativistic quantum chromodynamics (NRQCD). At leading
order in the velocity expansion the spin splittings are of $O(M_Q v^4)$, where
$M_Q$ is the renormalized quark mass and $v^2$ is the mean squared quark
velocity. A systematic analysis is done of all next-to-leading order
corrections. This includes the addition of $O(M_Q v^6)$ relativistic
interactions, and the removal of $O(a^2 M_Q v^4)$ discretization errors in the
leading-order interactions. Simulations are done for both S- and P-wave mesons,
with a variety of heavy quark actions and over a wide range of lattice
spacings. Two prescriptions for the tadpole improvement of the action are also
studied in detail: one using the measured value of the average plaquette, the
other using the mean link measured in Landau gauge. Next-to-leading order
interactions result in a very large reduction in the charmonium splittings,
down by about 60% from their values at leading order. There are further
indications that the velocity expansion may be poorly convergent for
charmonium. Prelimary results show a small correction to the hyperfine
splitting in the Upsilon system.
