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