I try to show that $$ I=\frac{\epsilon}{2\pi}\int^{\infty}_{-\infty}\mathrm{d}\varphi\int^{2\pi}_{0}\mathrm{d}\theta(e^{i\theta}+e^{-\varphi})e^{-2\epsilon\cosh(\varphi)+2\epsilon\cos(\theta)}=1 $$ for all $\epsilon>0$.
My idea was to express the integral with Bessel functions. This yields $$ I=-\pi\epsilon\big({J_1}(-2i\epsilon)H_0^{(1)}(2i\epsilon)+J_0(-2i\epsilon)H_1^{(1)}(2i\epsilon)\big) $$ where the $J'$s are Bessel functions of the first kind and $H$ are Hankel functions. However I don't see how to proceed from here. Are there any properties of the Bessel functions which might help? Or is there another way to see that $I=1$? Numerically the result of this integral is indeed $1$ independent of $\epsilon$.
Look at the Wronskian formula, Digital Library of Mathematical functions, 10.5.3. You'll need to get the arguments of the Bessel functions to be the same as that of the Hankel functions, and that will introduce a minus sign.
https://dlmf.nist.gov/10.5