I have the integral $$I=\int^{\infty}_0\frac{x^n}{e^x-1}\text{ d}x$$ for real $n>0$ that I want to evaluate only with contour integration.
(I already know the identity that $\Gamma(s)\zeta(s)=\int^{\infty}_0\frac{x^{s-1}}{\exp(x)-1}\text{ d}x$)
I let $$f(z)=\frac{z^n}{e^z-1}$$ noting that there is a removable singularity at $0$, and poles every $z=2\pi im$ for $m\in\mathbb{Z}, m\neq0$. As a result, I set up the following contour $\mathcal{C}$ shown below, where I would take the limit as $R\rightarrow+\infty$ and $\epsilon\rightarrow+0$.
By the Residue theorem, $$\oint_{\mathcal{C}}f(z)\text{ d}z=\int_{\text{Right}}+\int_{\text{Top}}+\int_{\text{Left}}+\int_{\text{Bottom}}+\int_rf(z)\text{ d}z=0$$
Parameterizing the integrals yields $$\int^{2\pi}_0\frac{(R+iy)^n}{e^R e^{iy}-1}i\text{ d}y + \int^{\epsilon}_{R}\frac{(2\pi i+x)^n}{e^{2\pi i}e^x-1}\text{ d}x + \int^{0}_{2\pi-\epsilon}\frac{(iy)^n}{e^{iy}-1}i\text{ d}y + \int^{R}_0\frac{x^n}{e^x-1}\text{ d}x + \int^{-\frac\pi 2}_{0}\frac{(2\pi i+\epsilon e^{i\theta})^n}{e^{2\pi i}e^{\epsilon e^{i\theta}}-1}i\epsilon e^{i\theta}\text{ d}\theta=0$$
Since $f(z)$ approaches $0$ as $\Re(z)$ becomes large, $\int_{\text{Right}}f(z)\text{ d}z=\int^{2\pi}_0\frac{(R+iy)^n}{e^r e^{iy}-1}i\text{ d}y$ should go to $0$. More rigorously, I moved the limit as $R\rightarrow +\infty$ inside (the function should uniformly converge I think so this is legal?) and indeed it would go to $0$.
$\int_{\text{Bottom}}f(z)\text{ d}z=I$ so we can ignore that for now. I'm thinking that we possibly could use a Laurent series for $\int_rf(z)\text{ d}z=\int^{-\frac\pi 2}_{0}\frac{(2\pi i+\epsilon e^{i\theta})^n}{e^{2\pi i}e^{\epsilon e^{i\theta}}-1}i\epsilon e^{i\theta}\text{ d}\theta$ since we're taking the limit $\epsilon\rightarrow +0$ but honestly I'm not sure how. If I directly take the limit $\epsilon\rightarrow +0$ of the integrand, Wolfram Alpha tells me it's equal to $(2\pi i)^n$, which would mean the entire integral is equal to $2^{n-1}(i\pi)^{n+1}$, but since it's such a complicated function I'm not sure if it uniformly converges on those bounds and thus if it is legal or not...
I also tried using the Cauchy-Schwartz Inequality on $\int_{\text{Left}}f(z)\text{ d}z=\int^{0}_{2\pi-\epsilon}\frac{(iy)^n}{e^{iy}-1}i\text{ d}y$ where $$\left|\int^{0}_{2\pi}\frac{(iy)^n}{e^{iy}-1}i\text{ d}y\right|\le \int^{0}_{2\pi-\epsilon}\frac{\left|i^n\right|\left|y^n\right|}{\left|e^{iy}-1\right|}|i|\text{ d}y = \int^{0}_{2\pi-\epsilon}\frac{y^n}{\left|e^{iy}-1\right|}\text{ d}y\le \int^{0}_{2\pi-\epsilon}\frac{y^n}{\left|\left|e^{iy}\right|-\left|1\right|\right|}\text{ d}y$$ where the last part follows from $\frac{1}{|a-b|}\le\frac{1}{||a|-|b||}$ but this makes the denominator $0$ which is illegal so I probably did something wrong here.
So in the end, I would have the following equation: $$I-\int^{\infty}_{0}\frac{(2\pi i+x)^n}{e^x-1}\text{ d}x - \int^{2\pi}_{0}\frac{(iy)^n}{e^{iy}-1}i\text{ d}y=2^{n-1} (i\pi)^{n+1}$$
And... I'm stuck. If anyone can help it would be really appreciated!
Thanks in advance! :D