  We study the spin splitting in the heavy quarkonium hybrid spectrum within
the framework of an nonrelativistic effective field theory. We derive for the
first time the spin-dependent part of the heavy-quark-antiquark potential for
heavy quarkonium hybrids to order $1/m^2$ in the heavy-quark-mass expansion. We
find that several operators that are not found in standard quarkonia appear,
most remarkably an operator suppressed by only one power of the heavy-quark
mass. By matching the weakly-coupled pNRQCD to the effective field theory in
the regime of short heavy-quark-antiquark distances, we work out the matching
coefficients of the spin-dependent operators, which are factorized into a
perturbative and a nonperturbative part. The nonperturbative part can be
expressed in terms of purely gluonic correlators. We fit the nonperturbative
parts of the matching coefficients to lattice data of the charmonium hybrid
spectrum and obtain results that respect the power counting. Using the obtained
nonperturbative pieces, we compute the bottomonium hybrid spectrum with the
spin-dependent potential, for which results from the lattice are still sparse.
