Evaluate the definite integral $\int^{\infty }_{0}\frac{x \,dx}{e^{x} -1}$ using contour integration

I would have guessed this was a duplicate, but I wasn't able to find another instance of this question during a cursory search.

Hint The denominator has period $2 \pi i$, which suggests using the following contour $\Gamma_{\epsilon, R}$, $0 < \epsilon < \pi$, $\epsilon < R$, (for which an illustration was already drawn for an answer to the similar question linked by Zacky in the comments):

The key trick here, which we apply with the benefit of hindsight, is to evaluate instead the similar integral $$\int_{\Gamma_{\epsilon, R}} \frac{z^2 \,dz}{e^z - 1} .$$ The interior of $\Gamma_{\epsilon, R}$ contains no poles, so this integral vanishes. Thus, parameterizing the constituent arcs of the contour gives \begin{multline} 0 = \underbrace{\int_\epsilon^R \frac{x^2 \,dx}{e^x - 1}}_{A} + \underbrace{\int_0^{2 \pi} \frac{(R + i y)^2 \cdot i \,dy}{e^{R + i y} - 1}}_{B} + \underbrace{\int_R^\epsilon \frac{(x + 2 \pi i)^2 \,dx}{e^x - 1}}_{C} \\ + \underbrace{\int_0^{-\pi / 2} \frac{(2 \pi i + \epsilon e^{i\theta})^2 \cdot i \epsilon e^{i \theta} d \theta}{e^{\epsilon e^{i \theta}} - 1}}_{D} + \underbrace{\int_{2 \pi - \epsilon}^\epsilon \frac{(i y)^2 \cdot i \,dy}{e^{i y} - 1}}_{E} + \underbrace{\int_{\pi / 2}^0 \frac{(\epsilon e^{i\theta})^2 \cdot i \epsilon e^{i \theta} d \theta}{e^{\epsilon e^{i \theta}} - 1}}_{F} . \qquad (\ast) \end{multline}

A standard bounding argument shows that $B \to 0$ as $R \to \infty$. Computing the first terms of the Taylor series gives that the integrand of $D$ is $-4 \pi^2 i + O(\epsilon)$, so $D = 2 \pi^3 i + O(\epsilon)$, and similarly $F = O(\epsilon)$ (in fact, the integrand is analytic at $0$, which implies this without any more computation). Now, expanding the integrand of $C$ gives $$-\int_\epsilon^R \frac{x^2 \,dx}{e^x - 1} = -\int_\epsilon^R \frac{x^2 \,dx}{e^x - 1} - 4 \pi i \int_\epsilon^R \frac{x \,dx}{e^x - 1} + 4 \pi^2 \int_\epsilon^R \frac{\,dx}{e^x - 1} .$$ The first term on the r.h.s. cancels $A$, and after taking appropriate limits the second term will be constant multiple of the integral $\color{#df0000}{\int_0^\infty \frac{x \,dx}{e^x - 1}}$ of interest. The third term diverges as $\epsilon \searrow 0$, and it turns out that the diverging part of this term in $\epsilon$ is canceled by the diverging part of $E$, but we can avoid dealing with this issue directly by passing to the imaginary part of $(\ast)$. Computing gives $\operatorname{Im} E = -\frac{1}{2} \int_\epsilon^{2 \pi - \epsilon} y^2 \,dy = -\frac{4}{3} \pi^3 + O(\epsilon)$, so taking the limits $\epsilon \searrow 0, R \to \infty$ of the imaginary part of $(\ast)$ leaves $$0 = -4 \pi \color{#df0000}{\int_0^\infty \frac{x\,dx}{e^x - 1}} + 2 \pi^3 - \frac{4}{3} \pi^3 ,$$ and rearranging gives the desired result, $$\color{#df0000}{\boxed{\int_0^\infty \frac{x \,dx}{e^x - 1} = \frac{\pi^2}{6}}} .$$

With contour integrals I think the easiest way is

$$\lim_{N \to \infty} 2\int_0^\infty x\frac{(1-e^{-Nx})}{e^x-1}dx =\lim_{N \to \infty} 2\sum_{n=1}^N \int_0^\infty xe^{-nx}dx=\lim_{N \to \infty}2\sum_{n=1}^N n^{-2}\\=\lim_{N \to \infty} \int_{|z| = N+1/2} \frac{z^{-2}}{e^{2i \pi z}-1}dz - 2i \pi Res(\frac{2i \pi z^{-2}}{e^{2i \pi z}-1},z=0)= -2i \pi Res(\frac{z^{-2}}{e^{2i\pi z}-1},z=0)$$

Travis has the right idea with the contour, but not quite the right auxiliary function. Instead, we use the contour $\gamma = ([0,R)\times\{0\})\cup(\{0\}\times[0,2\pi])\cup([0,R)\times\{2\pi\})\cup(\{R\}\times[0,2\pi])$ and integrate the function $f(z) = z(z-2\pi i)/(e^z-1)$ over this contour. The key is that $f(z)$ has no poles anywhere on the contour, so we don't need to worry about avoiding them to keep the integral well-defined. Since $f(z)$ also has no poles inside the contour, the integral around $\gamma$ is equal to zero. Thus, \begin{multline} \oint_\gamma f(z)dz = -\int_0^R \frac{x(x-2\pi i)}{e^x-1}dx + \int_0^{2\pi} \frac{(it)(it+2\pi i)}{e^{it}-1}idt \\ + \int_0^R \frac{(x+2\pi i)x}{e^{x+2\pi i}-1}dx - \int_0^{2\pi}\frac{(R+it)(R+it+2\pi i)}{e^{R+it}-1}idt = 0. \end{multline} It's not hard to see that the asymptotic behavior of the fourth integral is dominated by the $e^R$ in the denominator, so it will go to zero as $R\rightarrow\infty$. Noting that $e^{x+2\pi i} = e^x$, we take the that limit and do some simplification to get $$ 4\pi i\int_0^\infty \frac{x\,dx}{e^x-1} -i\int_0^{2\pi}\frac{t(t-2\pi)}{e^{it}-1}dt = 0. \Longrightarrow \int_0^\infty \frac{x\,dx}{e^x-1} = \frac{1}{4\pi}\int_0^{2\pi}\frac{t(t-2\pi)}{e^{it}-1}dt $$ Using $(e^{it}-1)^{-1} = -[1+i\,\mathrm{cot}(t/2)]/2$ and substituting $u = t - \pi$ in the imaginary integral gives $$\int_0^\infty \frac{x\,dx}{e^x-1} = \frac{1}{8\pi}\left[\int_0^{2\pi}t(2\pi-t)dt+i\int_{-\pi}^\pi (\pi^2-u^2)\tan\left(\frac{u}{2}\right)du\right]. $$ The $u$ integrand is clearly odd, so the integral is zero (as it has to be since the LHS is purely real). Since $\int_0^{2\pi}t(2\pi-t)dt = 4\pi^3/3$, we have $$ \int_0^\infty \frac{x\,dx}{e^x-1} = \frac{\pi^2}{6}. $$