Find the value of $\int_{0}^{\infty}\frac{x^3}{e^x-1}\ln(e^x - 1)\,dx$ I'm trying to figure out how to evaluate the following:
$$
J=\int_{0}^{\infty}\frac{x^3}{e^x-1}\ln(e^x - 1)\,dx
$$
I'm tried considering $I(s) = \int_{0}^{\infty}\frac{x^3}{(e^x-1)^s}\,dx\implies J=-I'(1)$, but I couldn't figure out what $I(s)$ was. My other idea was contour integration, but I'm not sure how to deal with the logarithm. Mathematica says that $J\approx24.307$. 
I've asked a similar question and the answer involved $\zeta(s)$ so I suspect that this one will as well. 
 A: Mathematica says that the answer is
$$\pi^2\zeta(3)+12\zeta(5)$$
I will try to figure out how this can be proven.

Added: Let me compute the 2nd integral in Ron Gordon's answer:
\begin{align}\int_{0}^{\infty}\frac{x^3 e^{-x}}{1-e^{-x}}\ln(1-e^{-x})\,dx
&=-\frac32\int_0^{\infty}x^2\ln^2(1-e^{-x})\,dx=\\&=-\frac32\left[\frac{\partial^2}{\partial s^2}\int_0^{\infty}e^{-sx}\ln^2(1-e^{-x})\,dx\right]_{s=0}=\\
&=-\frac32\left[\frac{\partial^2}{\partial s^2}\int_0^{1}t^{s-1}\ln^2(1-t)\,dt\right]_{s=0}=\\
&=-\frac32\left[\frac{\partial^4}{\partial s^2\partial u^2}\int_0^{1}t^{s-1}(1-t)^u\,dt\right]_{s=0,u=0}=\\
&=-\frac32\left[\frac{\partial^4}{\partial s^2\partial u^2}\frac{\Gamma(s)\Gamma(1+u)}{\Gamma(1+s+u)}\right]_{s=0,u=0}=\\
&=-\frac{1}{2}\left(\pi^2\psi^{(2)}(1)-\psi^{(4)}(1)\right).
\end{align}
To obtain the last expression, one should expand the ratio of gamma functions to 2nd order in $u$, then to expand the corresponding coefficient to 2nd order in $s$.
Then we can use that $\psi^{(2)}(1)=-2\zeta(3)$ and $\psi^{(4)}(1)=-24\zeta(5)$ (cf formula (15) here) to obtain the quoted result.
A: Using the change of variables $ u=e^{-x} $, we have 

$$\int_{0}^{\infty}\frac{x^3}{e^x-1}\ln(e^x - 1)\,dx  = \int _{0}^{1}\!{\frac { \left( \ln  \left( u \right)  \right) ^{3}
  \ln  \left( 1-u \right)  }{u-1
}}{du}- \int _{0}^{1}\!{\frac { \left( \ln  \left( u \right)  \right)^{4}  }{u-1
}}{du}. $$

Now, just apply the technique which has been used to find the exact solution in this problem and the result will follow.
A: How about pulling factors of $e^{-x}$ from both the denominator and log terms?  Then you end up with two separate integrals:
$$\int_0^{\infty}dx \frac{x^4 \, e^{-x}}{1-e^{-x}} + \int_0^{\infty}dx \frac{x^3 \, e^{-x}}{1-e^{-x}} \log{(1-e^{-x})}$$
In both cases, you Taylor expand the denominator in $e^{-x}$.  For the first integral, this results in
$$\sum_{k=0}^{\infty} \int_0^{\infty}dx\, x^4 \, e^{-(k+1) x} = 4! \sum_{k=0}^{\infty} \frac{1}{(k+1)^5} = 24 \, \zeta(5) $$
For the second integral, you also need to Taylor expand the log term.  This results in a double sum:
$$\begin{align}\sum_{k=0}^{\infty} \int_0^{\infty}dx\, x^3 \, e^{-(k+1) x} \log{(1-e^{-x})} &= -\sum_{k=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{m} \int_0^{\infty} dx \, x^3 e^{-(k+m) x}\\ &= - 3! \sum_{m=1}^{\infty} \frac{1}{m} \sum_{k=1}^{\infty} \frac{1}{(k+m)^4}\\ &= -\sum_{m=1}^{\infty} \frac{\psi^{(3)}(m+1)}{m} \end{align}$$
where $\psi$ is a polygamma function.
A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\on}[1]{\operatorname{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
J & \equiv \bbox[5px,#ffd]{\int_{0}^{\infty}{x^{3} \over
\expo{x} - 1}\ln\pars{\expo{x} - 1}\,\dd x}
\\[5mm] & =
\left.\partiald{}{\nu}\int_{0}^{\infty}x^{3}
\pars{\expo{x} - 1}^{\nu}\,\dd x
\,\right\vert_{\,\nu\ =\ -1}
\\[5mm] & =
\left.\partiald{}{\nu}\int_{0}^{\infty}x^{3}\expo{\nu x}
\pars{1 - \expo{-x}}^{\nu}\,\dd x
\,\right\vert_{\,\nu\ =\ -1}
\\[5mm] & =
\left.\partiald{}{\nu}\sum_{k = 0}^{\infty}{\nu \choose k}
\pars{-1}^{k}\int_{0}^{\infty}x^{3}\expo{-\pars{k - \nu}x}
\,\,\dd x\,\right\vert_{\,\nu\ =\ -1}
\\[5mm] & =
\left.\partiald{}{\nu}\sum_{k = 0}^{\infty}
{-\nu + k - 1\choose k}\,{6 \over \pars{k - \nu}^{4}}\,\right\vert_{\,\nu\ =\ -1}
\\[5mm] & =
24\sum_{k = 0}^{\infty}{1 \over \pars{k + 1}^{5}} -
6\sum_{ k = 0}^{\infty}{H_{k} \over \pars{k + 1}^{4}}
\\[5mm] & =
30\,\underbrace{\sum_{k = 1}^{\infty}{1 \over k^{5}}}
_{\ds{\zeta\pars{5}}}\ -\
6\
\underbrace{\sum_{ k = 1}^{\infty}{H_{k} \over k^{4}}}
_{\ds{3\zeta\pars{5} - \pi^{2}\zeta\pars{3}/6}}
\\[5mm] & =
\bbx{\pi^{2}\,\zeta\pars{3} + 12\,\zeta\pars{5}}
\approx 24.3070 \\ &
\end{align}
$\ds{\sum_{ k = 1}^{\infty}{H_{k} \over k^{4}}}$: See $\ds{\pars{20}}$ in MW.
