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

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}$.

Answer

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}}} .$$