A curious integral

The following integral has been on my mind for a while $$\int_0^\infty \frac{\sin(x)}{e^x-1}\,\mathrm d x \tag{$\dagger$}$$

Let us indicate the integrand as $f(x)=\frac{\sin(x)}{e^x-1}$. The following are a couple of observations.

1. The integrand can be extended by continuity in $0$ since $$\lim_{x\to 0}f(x)=\lim_{x\to 0}\frac{\frac{\sin(x)}{x}}{\frac{e^x-1}{x}}=1$$ This is the original reason I started playing around with this integral.

2. Mathematica yields the result $$\int_0^\infty \frac{\sin(x)}{e^x-1}\,\mathrm d x=\frac{\pi}{2}\textrm{Coth}(\pi)-\frac 12 \approx 1.076674047$$ which, following numerical evidence, seems correct.

3. Complex analysis may be useful here, since the integrand is a holomorphic function on $\mathbb C$. I tried writing $\sin(z)=\frac{e^{iz}-e^{-iz}}{2i}$, but I was not able to find an appropriate integration contour to solve the problem.

4. The $-1$ in the denominator breaks the simmetry of the expression. This made most of my substitutions useless.

Can this integral be evaluated correctly, preferrably through complex analytic methods?

Without resorting to the Abel-Plana formula, you may simply notice that

$$ I=\int_{0}^{+\infty}\frac{\sin x}{e^x-1}\,dx = \sum_{n\geq 1}\int_{0}^{+\infty}\sin(x)e^{-nx}\,dx = \sum_{n\geq 1}\frac{1}{n^2+1} \tag{1}$$ and since $\frac{\sin(\pi x)}{\pi x}=\prod_{n\geq 1}\left(1-\frac{x^2}{n^2}\right)$, by applying $\frac{d}{dx}\log(\cdot)$ to both sides we get: $$ -\frac{1}{x}+\pi\cot(\pi x) = \sum_{n\geq 1}\frac{2x}{x^2-n^2}\tag{2} $$ as well as: $$ \sum_{n\geq 1}\frac{1}{n^2+z^2} = \frac{-1+\pi z \coth(\pi z)}{2 z^2}\tag{3}$$ from which $I=\frac{-1+\pi\coth(\pi)}{2}$ clearly follows.
Have a look at this thread for further proofs.

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} \int_{0}^{\infty}{\sin\pars{x} \over \expo{x} - 1}\,\dd x & = {1 \over 2\ic}\int_{0}^{\infty}{\expo{-\ic x} - \expo{\ic x} \over 1 - \expo{-x}}\,\pars{-\expo{-x}}\,\dd x \,\,\,\stackrel{t\ =\ \expo{\large -x}}{=}\,\,\, -\,{1 \over 2\ic}\int_{0}^{1}{t^{\large\ic} - t^{\large -\ic} \over 1 - t} \,\dd t \\[5mm] & = -\,{1 \over 2\ic}\pars{H_{-\ic} - H_{\ic}}\qquad\qquad\qquad \pars{~H_{z}:\ Harmonic Number~} \\[5mm] & = -\,{1 \over 2\ic}\bracks{H_{-\ic} - \pars{H_{\ic - 1} + {1 \over \ic}}} \qquad\qquad\qquad\pars{~H_{z}\ Recurrence~} \\[5mm] & = -\,{1 \over 2} - {1 \over 2\ic}\bracks{H_{-\ic} - H_{\ic - 1}} = -\,{1 \over 2}-\,{1 \over 2\ic}\bracks{\pi\cot\pars{\pi\ic}} \pars{\substack{Euler \\[1mm] Reflection\ Formula}} \\[5mm] & = -\,{1 \over 2}-\,{\pi \over 2\ic}\bracks{-\ic\coth\pars{\pi}} = \bbx{{\pi \over 2}\,\coth\pars{\pi} - {1 \over 2}} \approx 1.077 \end{align}

Contour integration solution: We consider the integral

$$I = \oint_{C\left(\epsilon,R\right)} \frac{\exp(i z)dz}{\exp(z) - 1}$$

with $C\left(\epsilon,R\right)$ a rectangle from $\epsilon$ to $R$ on the real axis, from there to $R + 2\pi i$ parallel to the imaginary axis, from there to $\epsilon + 2\pi i$ parallel to the real axis, then a clockwise quarter turn with radius $\epsilon$ and center $2\pi i$ to the point $2\pi i - i\epsilon$, from there we move on the imaginary axis to the point $i\epsilon$ and then we take a clockwise quarter turn with radius $\epsilon$ and center the origin to move back to the starting point at $\epsilon$.

The integrand is analytic inside the contour, therefore $I = 0$. The sum of the two parts parallel to the real axis is:

$$I_r = \left[1-\exp(-2\pi)\right]\int_\epsilon^R \frac{\exp(i x)dx}{\exp(x) - 1}$$

So, the desired integral will follow from the imaginary part of $I_r$. The part of the contour integral from $R$ to $R + i$ tends to zero in the limit of $R\to\infty$ so this can be disregarded. The part of the contour integral along the imaginary axis can be written as:

$$I_i = -\int_{\epsilon}^{2\pi-\epsilon}\frac{\exp(-y)\exp\left(-i \frac{y}{2}\right)}{2\sin\left(\frac{y}{2}\right)}dy$$

We then see that

$$\lim_{\epsilon\to 0}\operatorname{Im}I_i = \frac{1-\exp(-2\pi)}{2}$$

The two quarter circles can be evaluated, we can borrow from the derivation of the residue theorem that each of them in the limit $\epsilon\to 0$ can be evaluated as $-\frac{\pi}{2} i$ times the residue at the poles that are at the centers of the quarter circles (unlike the case of a complete contour, this is only valid in the limit $\epsilon\to 0$). The residue at the pole at $z = 0$ is $1$, while the residue at $z = 2\pi i$ is $\exp(-2\pi)$. Having completed the evaluation of all parts of the contour integral, we can now equate the imaginary part of the contour integral to zero, and take the limit of $\epsilon\to 0$ and $R\to \infty$. This yields:

$$\int_0^{\infty}\frac{\sin(x)dx}{\exp(x)-1} = \frac{\pi}{2}\coth(\pi) - \frac{1}{2}$$

The Abel-Plana theorem states $$\sum _{n=0} ^{\infty} f(n) - \int_0^{\infty} f(s) \mathrm{d}s = \frac{1}{2} f(0) +i \int_0^{\infty} \frac{f(iy) - f(-iy)}{e^{2 \pi y} - 1} $$

Expressing the $sin$ using complex exponential your integral reads $$ \frac{1}{2i}\int_0^{\infty} \frac{e^{iz} - e^{-iz}}{e^z -1} \mathrm{d}z$$ by scaling the variable of integration, $$ \frac{2 \pi}{2i}\int_0^{\infty} \frac{e^{2 \pi iy} - e^{-2 \pi iy}}{e^{2 \pi y} -1} \mathrm{d}z = -\pi i\int_0^{\infty} \frac{e^{2 \pi iy} - e^{-2 \pi iy}}{e^{2 \pi y} -1} \mathrm{d}z$$

By stating $$ f(x) = - e^{-2 \pi x} $$ one can compute the series and the integral on the l.h.s of the Theorem

$$ \sum _{n=0} ^{\infty} f(n) - \int_0^{\infty} f(s) \mathrm{d}s = -\frac{e^{2 \pi i}}{ e^{2 \pi } -1} + \frac {1}{2 \pi} $$

Putting it all together, the sought integral equals

$$ - \pi [-\frac{e^{2 \pi i}}{ e^{2 \pi } -1} + \frac {1}{2 \pi} +\frac{1}{2}] = \frac{\pi}{2}[-1 +\frac{2 e^{2 \pi}}{e^{2 \pi}-1}]-\frac{1}{2}$$ which coincides with the given solution as $$ coth (x) = \frac{e^{2x} + 1}{e^{2x} - 1}$$