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

Question

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!

Answer 1

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