Help with logarithmic definite integral: $\int_0^1\frac{1}{x}\ln{(x)}\ln^3{(1-x)}\, dx$

I'm look for a closed form evaluation of the following improper definite integral involving logarithms:

$$\begin{align} I:&=\int_{0}^{1}\frac{1}{x}\ln{(x)}\ln^3{(1-x)}\,\mathrm{d}x\\ &\approx 0.579307. \end{align}$$

My attempt:

My first idea was to transform the integral using the substitution $u=-\ln{(x)}$, but this didn't really result in anything recognizably easier than the original integral:

$$\begin{align} I&=\int_{0}^{1}\frac{1}{x}\ln{(x)}\ln^3{(1-x)}\,\mathrm{d}x\\ &=\int_{0}^{\infty}(-u)\ln^3{(1-e^{-u})}\,\mathrm{d}u\\ &=-\int_{0}^{\infty}u\ln^3{(1-e^{-u})}\,\mathrm{d}u. \end{align}$$

The second thing I tried was a substitution I've known to be successful with similar problems:

$$\begin{align} u&=-\ln{(1-x)},\\ -u&=\ln{(1-x)},\\ e^{-u}&=1-x,\\ x&=1-e^{-u},\\ \mathrm{d}x&=e^{-u}\,\mathrm{d}u. \end{align}$$

Applying the substitution, the integral becomes:

$$\begin{align} I&=\int_{0}^{1}\frac{1}{x}\ln{(x)}\ln^3{(1-x)}\,\mathrm{d}x\\ &=\int_{0}^{\infty}\frac{(-u)^3}{1-e^{-u}}\ln{(1-e^{-u})}\,e^{-u}\,\mathrm{d}u\\ &=-\int_{0}^{\infty}\frac{u^3\,e^{-u}}{1-e^{-u}}\ln{(1-e^{-u})}\,\mathrm{d}u\\ &=\int_{0}^{\infty}\frac{u^3}{1-e^{u}}\ln{(1-e^{-u})}\,\mathrm{d}u. \end{align}$$

This looks like it has a shot at being evaluated in terms of series expansions, but I'm stuck at this point.

Suggestions welcome.

Update:

Using the series expansions $\ln{(1-e^{-x})}=-\sum_{n=1}^{\infty}\frac{1}{n}e^{-nx}$ and $\sum_{k=1}^{\infty}e^{-kx}=\frac{1}{e^{x}-1}$, the integral becomes:

$$\begin{align} I&=\int_{0}^{\infty}\frac{x^3}{1-e^{x}}\ln{(1-e^{-x})}\,\mathrm{d}x\\ &=\int_{0}^{\infty}\frac{x^3}{e^{x}-1}\sum_{n=1}^{\infty}\frac{1}{n}e^{-nx}\,\mathrm{d}x\\ &=\sum_{n=1}^{\infty}\frac{1}{n}\int_{0}^{\infty}\frac{x^3\,e^{-nx}}{e^{x}-1}\,\mathrm{d}x\\ &=\sum_{n=1}^{\infty}\frac{1}{n}\int_{0}^{\infty}x^3\,e^{-nx}\sum_{k=1}^{\infty}e^{-kx}\,\mathrm{d}x\\ &=\sum_{n=1}^{\infty}\frac{1}{n}\sum_{k=1}^{\infty}\int_{0}^{\infty}x^3\,e^{-(k+n)x}\,\mathrm{d}x\\ &=\sum_{n=1}^{\infty}\frac{1}{n}\sum_{k=1}^{\infty}\frac{6}{(k+n)^4}\\ &=\sum_{n=1}^{\infty}\frac{1}{n}\psi^{(3)}(n+1), \end{align}$$

where $\psi^{(n)}(x)$ is the $n$-th derivative of the digamma function. Unfortunately, I don't know how to evaluate this last series. But this feels like progress.

We could use the series expansion $$\frac{\ln(1-x)}{1-x}=-\sum_{n=1}^\infty H_nx^n,$$ where $H_n=1+\frac12+\frac13+\ldots+\frac1n$ is the harmonic number. Then we get $$I=\int_0^1\frac{\ln(1-x)\ln^3x}{1-x}dx=-\sum_{n=1}^\infty H_n\int_0^1 x^n\ln^3 x\,dx=6\sum_{n=1}^\infty \frac{H_n}{(n+1)^4}.$$ The last series could be evaluated using Euler's formula $$\sum_{n=1}^\infty \frac{H_n}{n^4}=3\zeta(5)-\zeta(2)\zeta(3)$$ (see formula (20) here: http://mathworld.wolfram.com/HarmonicNumber.html): $$I=6\sum_{n=1}^\infty \frac{H_n}{(n+1)^4}=6\sum_{n=1}^\infty \left(\frac{H_{n+1}}{(n+1)^4}-\frac{1}{(n+1)^5}\right)=6(2\zeta(5)-\zeta(2)\zeta(3)).$$

Introduce a two-parameter deformation $$\mathcal{I}(a,b)=\int_0^1 x^{a-1}\left[(1-x)^b-1\right] dx=\frac{\Gamma(a)\Gamma(b+1)}{\Gamma(a+b+1)}-\frac1a.$$ The integral we are looking for is obtained as \begin{align*} I&=\frac{\partial^3}{\partial b^3}\left[\frac{\partial\mathcal{I}(a,b)}{\partial a}\biggl|_{a=0}\right]_{b=0}=\\ &=\frac{\partial^3}{\partial b^3}\left[\frac{\gamma^2}{2}+\frac{\pi^2}{12}+\gamma\,\psi(1+b)+\frac{\psi^2(1+b)-\psi'(1+b)}{2}\right]_{b=0}=\\ &=12\zeta(5)-6\zeta(2)\zeta(3). \end{align*} At the last step, one needs to use that $\psi^{(n)}(1)=(-1)^{n+1}n!\,\zeta(n+1)$.

If you're interested in evaluating that $\psi'''(n+1)$ series you obtained in your original evaluation, a fun way to go about it is to use contours. A method popularized by Flajolet and Salvy in their famous paper on the topic.

$$\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}$$

Begin with the kernel $$f(z)=\pi \cot(\pi z)\psi'''(-z)$$

This has poles at the positive integers, n; negative integers, -n, and 0.

The series at the positive integers is:

$$\frac{6}{(z-n)^{5}}-\frac{2\pi^{2}}{(z-n)^{3}}-\frac{\psi'''(n+1)}{z-n}+\cdot\cdot\cdot $$

Thus, the sum of the residues is:

$$Res(f,n)=\frac{6}{4!}\frac{d^{4}}{dz^{4}}\left[\frac{1}{(z-n)^{5}}\cdot \frac{1}{n}\right]-\frac{2\pi^{2}}{2!}\frac{d^{2}}{dz^{2}}\left[\frac{1}{(z-n)^{3}}\cdot \frac{1}{n}\right]-\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}$$

$$=6\sum_{n=1}^{\infty}\frac{1}{n^{5}}-2\pi^{2}\sum_{n=1}^{\infty}\frac{1}{n^{3}}-\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}$$

$$=6\zeta(5)-2\pi^{2}\zeta(3)-\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}$$

the series at the negative integers has a simple pole.

$$Res(f, -n)=\sum_{n=1}^{\infty}\frac{\psi'''(n)}{n}$$$

The residue at z=0 can be found by the Laurent expansion and noting the coefficient of the 1/z term:

$$Res(f,0)=24\zeta(5)$$

Now, put it altogether and set equal to 0 due to there being no poles in the contour:

$$6\zeta(5)-2\pi^{2}\zeta(3)+24\zeta(5)+\sum_{n=1}^{\infty}\frac{\psi'''(n)}{n}-\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}=0$$

solve for the series in question. Note, since $\frac{\psi'''(n)-\psi'''(n+1)}{n}=\frac{6}{n^{5}}$, this means that $\displaystyle \sum_{n=1}^{\infty}\frac{\psi'''(n)}{n}$ is equal to the series in question plus $6\zeta(5)$. In other words by calling our series S:

$\displaystyle \sum_{n=1}^{\infty}\frac{\psi'''(n)}{n}=S+6\zeta(5)$

$$-2S+24\zeta(5)-2\pi^{2}\zeta(3)=0$$

$$S=12\zeta(5)-\pi^{2}\zeta(3)$$