calculation of Stefan's constant In the calculation of Stefan's constant one has the integral
$$J=\int_0^\infty \frac{x^{3}}{\exp\left(x\right)-1} \, dx$$
which according to Wikipedia is equal to $\frac{\pi^4}{15}$. 
In this page of Wikipedia there is a (long) method of calculation using  the Taylor expansion of $f(k) = \int_0^\infty \frac{\sin\left(kx\right)}{\exp\left(x\right)-1} \, dx$ obtained with contour integration.
In this article is written that $J=\Gamma(4)\,\mathrm{Li}_4(1) = 6\,\mathrm{Li}_4(1) = 6 \zeta(4)$. How to obtain these equalities?
It is also written that there is number way to obtain the result. Someone know another way?
 A: Writing the integral as 
$$
   J = \int_0^\infty \frac{x^3}{\mathrm{e}^{x} -1} \mathrm{d}x = \int_0^\infty \frac{x^3 \mathrm{e}^{-x} }{ 1 - \mathrm{e}^{-x}} \mathrm{d}x = \lim_{\epsilon \downarrow 0} \int_\epsilon^\infty \frac{x^3 \mathrm{e}^{-x} }{ 1 - \mathrm{e}^{-x}} \mathrm{d}x
$$
Now for $x >\epsilon$, $\left(1-\mathrm{e}^{-x}\right)^{-1} = \sum_{k=0}^\infty \mathrm{e}^{-k x}$. Interchanging the summation and integration, justified by Tonelli's theorem:
$$ \begin{eqnarray}
  J &=& \lim_{\epsilon \downarrow 0} \sum_{k=0}^\infty \int_{\epsilon}^\infty x^3 \mathrm{e}^{-(k+1)x} \mathrm{d}x \\  &=& \lim_{\epsilon \downarrow 0} \sum_{k=0}^\infty \frac{\Gamma(4)}{(k+1)^4} \mathrm{e}^{-k \epsilon} \left(1 + (k+1) \epsilon + \frac{1}{2} (k+1)^2 \epsilon^2 + \frac{1}{3!} (k+1)^3 \epsilon^3 \right) \\ &=& \Gamma(4) \lim_{\epsilon \downarrow 0} \left(  \operatorname{Li}_4(\mathrm{e}^{-\epsilon}) + \epsilon \operatorname{Li}_3(\mathrm{e}^{-\epsilon}) + \frac{\epsilon^2}{2!} \operatorname{Li}_2(\mathrm{e}^{-\epsilon}) + \frac{\epsilon^3}{3!} \operatorname{Li}_1(\mathrm{e}^{-\epsilon}) \right)
\end{eqnarray}
$$
Since $\operatorname{Li}_1(\mathrm{e}^{-\epsilon}) = -\log\left(1-\mathrm{e}^{-\epsilon}\right)$ and $\lim_{\epsilon \downarrow 0} \epsilon^3 \log\left(1-\mathrm{e}^{-\epsilon}\right) = 0$, and due to finiteness of $\operatorname{Li}_{k}(1)$ for $k \geqslant 2$ we have
$$
   J = \Gamma(4) \operatorname{Li}_4(1) = \Gamma(4) \zeta(4) = \frac{\pi^4}{15}
$$
A: Let us compute something more general the integrals defined by
$$ I_{n,p}=\int_0^\infty \frac{x^n}{(1-\mathrm{e}^{-x})^{p+1}}\mathrm dx$$
that appear in statistical physics when studying bosonic systems.
The simplest method is to develop the denominator into a series
$$ \frac{1}{(1-\mathrm e^{-x})^{p+1}}=
\sum_{k=0}^\infty \binom{k+p}p\left(\mathrm e^{-x}\right)^k $$
and use the integral 
$$ \int_0^\infty x^n\mathrm{e}^{-kx}\mathrm dx= \frac{n!}{k^{n+1}}.$$
Combining these, one gets
$$I_{n,p}=\frac{n!}{p!} \sum_{k=0}^\infty\frac{(k+p)(k+p-1)\cdots(k+1)}{k^{n+1}}\tag{1}$$
in which we recognize the definition of Riemann's zeta function
$\zeta(k)=\sum_{n=1}^\infty n^{-k}$. The numerator in (1) can be rewritten using Stirling numbers of the first kind thanks to
$$ (k+p)\cdots(k+1)=\sum_{m=0}^{p}\genfrac{[}{]}{0pt}{}{p+1}{m+1}k^m$$
and we obtain finally a closed form
$$ I_{n,p}=\frac{n!}{p!}\sum_{m=0}^p\genfrac{[}{]}{0pt}{}{p+1}{m+1}\zeta(n+1-m).$$
In the case $n=3$, $p=0$ the result is
$$I_{3,0}=3!\,\times\zeta(4)=6\times\frac{\pi^4}{90}=\frac{\pi^4}{15}.$$
