Is the eigenvalue decomposition of the Sturm-Liouville operator $$ Lf(x)=-f''(x)+h\sin(x)f'(x),\quad h>0, $$ with Neumann boundary conditions $f'(-\pi)=f'(\pi)=0$ on the Hilbert space $L^2([-\pi,\pi],\mu_h)$ known? Here $d\mu_h(x)=Z^{-1}e^{h\cos x}dx$ and $Z>0$ is chosen such that $\mu_h$ is a probability measure. I suspect it involves Bessel functions.

N.B. The Sturm-Liouville problem is unitarily equivalent to the Schrödinger operator $$ Hf(x)=-f''(x)+\frac{h}{4}(h\sin^2(x)-2\cos(x))f(x) $$ on $L^2([-\pi,\pi])$.

**Edit:** By the hint of Sascha, we transform $H$ into a Schrödinger operator with Whittaker-Hill potential. Set $y=x/2$, then
$$
Hf(y)=-\frac14 f''(y)+\left(\frac{h^2}{8}-\frac{h^2}{8}\cos(4y)-\frac{h}{2}\cos(2y)\right)f(y).
$$
Hence, $4H-h^2/2$ is a Schrödinger operator with Whittaker-Hill potential with parameters $\alpha=h/2$, $s=1$ (in the convention of the paper by Hemery and Veselov).

Is anything known about the spectral gap of this operator (as noted by Hemery and Veselov, the ground state eigenvalue is (not very surprisingly) $-h^2/2$)?