How can I find the formula for the integral $\int_{0}^{\infty} \frac{x^{n}\left(e^{3 x}-e^{x}\right)}{\left(e^{x}-1\right)^{4}} d x ?$ where $n>2$.

Couple of days ago, I was asked to find the exact value of the definite integral

$$\displaystyle \int_{0}^{\infty} \frac{x^{5}\left(e^{3 x}-e^{x}\right)}{\left(e^{x}-1\right)^{4}} d x ,\tag*{}$$ After finding its exact value, I noted that the proof can be generalized to $$n$$ instead of 5. Now I am going to share my proof with you all.

$$\displaystyle I(n)=\int_{0}^{\infty} \frac{x^{n}\left(e^{3 x}-e^{x}\right)}{\left(e^{x}-1\right)^{4}} d x,$$ where $$n\geq 3.$$

Simplifying and splitting the integrand into 2 simpler one yields $$\displaystyle I(n)=\int_{0}^{\infty} \frac{x^{n} e^{x}\left(e^{x}+1\right)}{\left(e^{x}-1\right)^{3}} d x=\underbrace{\int_{0}^{\infty} \frac{x^{n} e^{x}}{\left(e^{x}-1\right)^{2}} d x}_{J} +2\underbrace{\int_{0}^{\infty} \frac{x^{n} e^{x}}{\left(e^{x}-1\right)^{3}} d x}_{K} \tag*{}$$ Using integration by parts and geometric series, we evaluate the integral $$J.$$ \displaystyle \begin{aligned}J &=-\int_{0}^{\infty} x^{n} d\left(\frac{1}{e^{x}-1}\right) \\&=-\left[\frac{x^{n}}{e^{x}-1}\right]_{0}^{\infty}+n\int_{0}^{\infty} \frac{x^{n-1}}{e^{x}-1} d x \\&=n\int_{0}^{\infty} \frac{x^{n-1} e^{-x}}{1-e^{-x}} d x \\&=n\int_{0}^{\infty} x^{n-1} e^{-x}\left(\sum_{k=0}^{\infty} e^{-k x}\right) d x\\&=n\sum_{k=0}^{\infty} \int_{0}^{\infty} x^{n-1} e^{-(k+1) x} d x\end{aligned} \tag*{} Similarly for the integral $$K$$, we have \displaystyle \begin{aligned}K &=\int_{0}^{\infty} \frac{x^{n} e^{x}}{\left(e^{x}-1\right)^{3}} d x \\&=-\frac{1}{2} \int_{0}^{\infty} x^{n} d\left(\frac{1}{e^{x}-1}\right)^{2} \\&=\frac{n}{2} \int_{0}^{\infty} \frac{x^{n-1}}{\left(e^{x}-1\right)^{2}} d x \\&=\frac{n}{2} \int_{0}^{\infty} \frac{x^{n-1} e^{-2 x}}{\left(1-e^{-x}\right)^{2}} d x\end{aligned} \tag*{} To deal with the last integral, we do need an infinite geometric series $$\displaystyle \frac{1}{1-y}=\sum_{k=0}^{\infty} y^{k} \text { for }|y|<1.\tag*{}$$ Now we differentiate both sides w.r.t. $$y$$ and obtain $$\displaystyle \frac{1}{(1-y)^{2}}=\sum_{k=0}^{\infty} k y^{k-1} \tag*{}$$ Replacing $$y$$ by $$e^{-x}$$yields \begin{aligned}K &=\frac{n}{2} \int_{0}^{\infty} x^{n-1} e^{-2 x} \sum_{k=0}^{\infty} k e^{-(k-1) x} d x \\&=\frac{n}{2} \sum_{k=0}^{\infty} k\int_{0}^{\infty} x^{n-1} e^{-(k+1) x} d x\end{aligned}\tag*{} Grouping them together, we can conclude that \displaystyle \begin{aligned}I(n) &=J+2 K \\&=n\left(\sum_{k=0}^{\infty} \int_{0}^{\infty} x^{n-1} e^{-(k+1) x} d x+\sum_{k=0}^{\infty} k \int_{0}^{\infty} x^{k} e^{-(k+1) x} d x\right)\\&=n \sum_{k=0}^{\infty}(k+1) \int_{0}^{\infty} x^{n-1} e^{-(k+1) x} d x \\&=n\sum_{k=0}^{\infty}(k+1) \cdot \frac{(n-1)!}{(k+1)^{n}} \quad \text {Via IBP repeatedly}\\&=n! \sum_{k=1}^{\infty} \frac{1}{k^{n-1}}\end{aligned} \tag*{} I finally succeed to find a beautiful formula for the integral $$I(n)$$ with $$n\geq 3$$, $$\boxed{\int_{0}^{\infty} \frac{x^{n}\left(e^{3 x}-e^{x}\right)}{\left(e^{x}-1\right)^{4}} d x =n! \zeta (n-1) }$$ :|D Wish you enjoy my proof! Your suggestions, comments and alternate methods are warmly welcome!

• I enjoy ! Well done and $\to +1$ Nov 6, 2021 at 9:02
• Solution was elegant $(+1)$ Nov 6, 2021 at 9:04
• Happy to hear from all of you. Thank you for your encouragement.
– Lai
Nov 6, 2021 at 9:11
• (+1) The general answer for any $b\in R; b>2$: $$\frac{d^2}{dx^2}\frac{1}{e^x-1}=\frac{2e^{2x}}{(e^x-1)^3}-\frac{e^x}{(e^x-1)^2}=\frac{e^{2x}+e^x}{(e^x-1)^3}=\frac{e^{3x}-e^x}{(e^x-1)^4}$$ Therefore, $$I(n)=\int_{0}^{\infty} \frac{x^{b}\left(e^{3 x}-e^{x}\right)}{\left(e^{x}-1\right)^{4}} d x=\displaystyle\int_0^\infty x^b\frac{d^2}{dx^2}\Big(\frac{1}{e^x-1}\Big)dx$$ Integrating by part two times $$=b(b-1)\int_0^\infty \frac{x^{b-2}}{e^x-1}dx =\Gamma(b+1)\zeta(b-1)$$ Nov 6, 2021 at 9:18
• You are right! I am sorry. Your solution is very elegant. Thank you for your generalization.
– Lai
Nov 6, 2021 at 9:29

$$I(n, p)= \int_0^{\infty}\frac {x^ne^x}{(e^x-1)^p}dx = \frac {n\Gamma(n)}{p-1}\sum_{k =0}\binom k{p-2}\frac 1{(k+1)^n}$$ $$I(n, 2) + 2\times I(n, 3) = n\Gamma(n)\zeta(1) + n\Gamma(n)(\zeta(n-1) -\zeta(1)) =n\Gamma(n)\times \zeta(n-1)$$
Let: $$I(5, 5) = -30\zeta(5) +55\zeta(4)-30\zeta(3)+5\zeta(2) + \color{red}{0}\times\zeta(1)$$ $$I(6, 6) = 144\zeta(6) - 300\zeta(5) +210\zeta(4)-60\zeta(3) + 6\zeta(2) + \color{red}{0}\times\zeta(1)$$
If integration converges: $$C_0\times\zeta(1)\implies$$ $$C_0 = 0$$