How to do integral $\int_{0}^{\infty} \frac{x^4 e^{x}}{(e^x-1)^2}\,dx$ $$\int_{0}^{\infty} \frac{x^4 e^{x}}{(e^x-1)^2}\,dx$$
This integral came from the Debye theory of the molecular heat capacity of crystal.
I heard that this integral form is related to Riemann Zeta function and Gamma function.
Also, i heard that it can be solved by using Feynmann's Trick.
How can i evaluate this integral?
 A: The Riemann/Gamma way:$$\int_0^\infty\frac{x^4e^{-x}dx}{(1-e^{-x})^2}=\sum_{n\ge 1}n\int_0^\infty x^4e^{-nx}dx=\sum_{n\ge 1}\frac{4!}{n^4}=24\zeta(4)=\frac{4\pi^4}{15}.$$
A: Recall the integral representation of the Polylogarithm given by

$$\operatorname{Li}_s(z)~=~\frac1{\Gamma(s)}\int_0^\infty \frac{t^{s-1}}{e^t/z-1}\mathrm dt$$

Moreover note that
$$\frac{\partial}{\partial z}\left[\frac{t^{s-1}}{e^t/z-1}\right]=\frac{t^{s-1}e^t}{(e^t-z)^2}$$
which is precisely the structure of your integrand. Thus, by applying a variation of Feynman's Trick we obtain
\begin{align*}
\frac{\mathrm d}{\mathrm dz}\Gamma(s)\operatorname{Li}_s(z)&=\frac{\mathrm d}{\mathrm dz}\int_0^\infty \frac{t^{s-1}}{e^t/z-1}\mathrm dt\\
\Gamma(s)\frac{\mathrm d}{\mathrm dz}\operatorname{Li}_s(z)&=\int_0^\infty \frac{\partial}{\partial z}\frac{t^{s-1}}{e^t/z-1}\mathrm dt\\
\Gamma(s)\frac{\operatorname{Li}_{s-1}(z)}z&=\int_0^\infty\frac{t^{s-1}e^t}{(e^t-z)^2}\mathrm dt
\end{align*}
Now plugging in $s=5$ and $z=1$ yields to
\begin{align*}
\Gamma(5)\frac{\operatorname{Li}_{5-1}(1)}1&=\int_0^\infty\frac{t^{5-1}e^t}{(e^t-1)^2}\mathrm dt\\
4!\operatorname{Li}_4(1)&=\int_0^\infty\frac{t^{4}e^t}{(e^t-1)^2}\mathrm dt
\end{align*}
Note that the RHS is the integral we are interested in whereas the LHS can be evaluated directly in terms of the Riemann Zeta Function. Overall this boils down to utilizing the value of $\zeta(4)$, as J.G.'s answer did aswell, since $\operatorname{Li}_s(1)=\zeta(s)$ as one can verify by using the series representation of the Polylogarithm.

$$\therefore~\int_0^\infty\frac{t^{4}e^t}{(e^t-1)^2}\mathrm dt~=~4!\operatorname{Li}_4(1)~=~\frac{4\pi^4}{15}$$

A: For a slight variation along a polylogarithmic theme, observe that for the polylogarithm of order $s = -1$ one has
$$\operatorname{Li}_{-1} (z) = \frac{z}{(1 - z)^2}.$$
Setting $z = e^{-x}$ leads to
$$\operatorname{Li}_{-1} (e^{-x}) = \frac{e^{-x}}{(1 - e^{-x})^2}.$$
Now for the integral
$$I = \int_0^\infty \frac{x^4 e^x}{(e^x - 1)^2} \, dx,$$
rewriting it as
$$I = \int_0^\infty \frac{x^4 e^{-x}}{(1 - e^{-x})^4} \, dx,$$
which in terms of polylogarithms becomes
$$I = \int_0^\infty x ^4 \operatorname{Li}_{-1} (e^{-x}) \, dx. \tag1$$
Now since
$$\int \operatorname{Li}_s (e^{-x}) \, dx = - \operatorname{Li}_{s + 1} (e^{-x}) + C,$$
integrating (1) four times by parts leads to
\begin{align}
I &= 4! \int_0^\infty \operatorname{Li}_3 (e^{-x}) \, dx = 4! \big{[} - \operatorname{Li}_4 (e^{-x}) \big{]}_0^\infty = 4! \operatorname{Li}_4 (1) = 4! \cdot \frac{\pi^4}{90},
\end{align}
or
$$\int_0^\infty \frac{x^4 e^x}{(e^x - 1)^2} \, dx = \frac{4 \pi^4}{15},$$
as required. 

Comment
The real advantage of this approach lay in the fact it can be used to find the corresponding indefinite integral 
$$\int \frac{x^4 e^x}{(e^x - 1)^2} \, dx.$$
Following exactly the method outlined above, one finds
$$\int \frac{x^4 e^{x}}{(e^x - 1)^2} \, dx = -x^4 \operatorname{Li}_0 (e^{-x}) - 4x^3 \operatorname{Li}_1 (e^{-x}) - 12 x^2 \operatorname{Li}_2 (e^{-x}) - 24x \operatorname{Li}_3 (e^{-x}) - 24 \operatorname{Li}_4 (e^{-x}) + C.$$
Note the polylogarithm expressions of order zero and one can be expressed as
$$\operatorname{Li}_0 (e^{-x}) = \frac{1}{e^x - 1} \qquad \text{and} \qquad \operatorname{Li}_1 (e^{-x}) = - \ln(1 - e^{-x}).$$ 
