  We construct a nonrelativistic effective field theory description of heavy
quarkonium hybrids from QCD. We identify the symmetries of the system made of a
heavy quark, a heavy antiquark, and glue in the static limit. Corrections to
this limit can be obtained order by order in an expansion in the inverse of the
mass $m$ of the heavy quark. At order $1/m$ in the expansion, we obtain at the
level of potential Non-Relativistic QCD a system of coupled Schr\"odinger
equations that describes hybrid spin-symmetry multiplets, including the mixing
of different static energies into the hybrid states, an effect known as
$\Lambda$-doubling in molecular physics. In the short distance, the static
potentials depend on two nonperturbative parameters, the gluelump mass and the
quadratic slope, which can be determined from lattice calculations. We adopt a
renormalon subtraction scheme for the calculation of the perturbative part of
the potential. We numerically solve the coupled Schr\"odinger equations and
obtain the masses for the lowest lying spin-symmetry multiplets for $c\bar{c}$,
$b\bar{c}$, and $b\bar{b}$ hybrids. The $\Lambda$-doubling effect breaks the
degeneracy between opposite parity spin-symmetry multiplets and lowers the mass
of the multiplets that get mixed contributions of different static energies. We
compare our findings to the experimental data, direct lattice computations, sum
rules calculations, and discuss the relation to the Born-Oppenheimer
approximation.
