My friend and I have been trying weeks to evaluate the integral
$$\int^{\infty }_{0}\frac{x \,dx}{e^{x} -1} .$$
We have together tried 23 contours, and all have failed.
We already know how to solve this with infinite sums (i.e., using the zeta function and the Basel problem), but we can't figure out how to solve it using contour integration methods.
We already know the answer is $\frac{\pi^{2}}{6}$.
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}. $$
