  We compute the heavy quarkonium mass of $l\not= 0$ (angular momentum) states,
with otherwise arbitrary quantum numbers, with next-next-to-next-to-leading
logarithmic (N$^3$LL) accuracy. This constitutes the first observable in heavy
quarkonium for which two orders of the weak-coupling expansion sensitive to the
ultrasoft scale are known and the resummation of ultrasoft logarithms is made.
We also obtain, for the first time, resummed N$^3$LL expressions for the
different fine and hyperfine energy splittings of these states, which are not
sensitive to the ultrasoft scale but still require resummation of (hard)
logarithms. We do this analysis for the equal and non-equal mass cases. We also
study an alternative computational scheme that treats the static potential
exactly. We then perform a comprehensive phenomenological analysis: we apply
these results to the $n=2$, $l=1$ bottomonium, $B_c$ and charmonium systems and
study their convergence.
