How to prove this series identity $\sum_{n=1}^{\infty}\left(\frac{1}{n}\sum_{k=n+1}^{\infty}\frac{1}{k^3}\right)=\frac{\pi^{4}}{360}$? I tried to prove this identity seemingly related to special functions which has been verified via Mathematica without a proof:
$$\sum_{n=1}^{\infty}\left(\frac{1}{n}\sum_{k=n+1}^{\infty}\frac{1}{k^3}\right)=\frac{\pi^{4}}{360}$$
I rewrote the left side as a improper integral:
$$\frac{1}{\Gamma(3)}\int_{0}^{\infty}\frac{x^2}{e^{x}-1}\ln(\frac{e^{x}}{e^{x}-1})\mathrm{d}x$$
But I have trouble in calculating
$$\int_{0}^{\infty}\frac{x^2}{e^{x}-1}\ln({e^{x}-1})\mathrm{d}x$$
Could anybody tell me how to calculate the last integral or find another approach to the original problem?
 A: By double counting 

we have
$$\sum_{n=1}^{\infty}\left(\frac{1}{n}\sum_{k=n+1}^{\infty}\frac{1}{k^3}\right)
=\sum_{n=2}^{\infty}\sum_{k=1}^{n-1}\frac{1}{kn^3}=$$
$$=\sum_{n=2}^{\infty}\left(\frac{1}{n^3}\sum_{k=1}^{n}\frac{1}{k}\right)-\sum_{n=2}^{\infty}\frac{1}{n^4}
=\sum_{n=2}^{\infty}\frac{H_{n}}{n^3}-\sum_{n=2}^{\infty}\frac{1}{n^4}=\frac{\pi^4}{72}-\frac{\pi^4}{90}=\frac{\pi^4}{360}$$
where we have used the result


*

*A series involving the harmonic numbers : $\sum_{n=1}^{\infty}\frac{H_n}{n^3}$
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{\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}
&\bbox[10px,#ffd]{\sum_{n = 1}^{\infty}\pars{{1 \over n}
\sum_{k = n + 1}^{\infty}{1 \over k^{3}}}} =
\sum_{k = 1}^{\infty}{1 \over k^{3}}\sum_{n = 1}^{k - 1}{1 \over n} =
\sum_{k = 1}^{\infty}{1 \over k^{3}}\pars{-\,{1 \over k} + \sum_{n = 1}^{k}{1 \over n}}
\\[5mm] = &
-\ \overbrace{\sum_{k = 1}^{\infty}{1 \over k^{4}}}
^{\zeta\pars{4}\ =\ \pi^{\large 4}/90}\ +\ 
\sum_{k = 1}^{\infty}{1 \over k^{3}}
\sum_{n = 1}^{\infty}\pars{{1 \over n} - {1 \over n + k}} =
-\,{\pi^{4} \over 90} + 
\sum_{k = 1}^{\infty}\sum_{n = 1}^{\infty}
{1 \over k^{2}n\pars{n + k}}
\\[5mm] = &\
-\,{\pi^{4} \over 90} + 
{1 \over 2}\sum_{k = 1}^{\infty}\sum_{n = 1}^{\infty}
\bracks{{1 \over k^{2}n\pars{n + k}} + {1 \over n^{2}k\pars{k + n}}} =
-\,{\pi^{4} \over 90} + 
{1 \over 2}\sum_{k = 1}^{\infty}\sum_{n = 1}^{\infty}{1 \over k^{2}n^{2}}
\\[5mm] = &\
-\,{\pi^{4} \over 90} + 
{1 \over 2}\pars{\pi^{2} \over 6}^{2} =
\bbx{\pi^{4} \over 360} \approx 0.2706 \\ &
\end{align}
A: A continuation to your work: 
$$\int_0^\infty\frac{x^2}{e^x-1}\ln\left(\frac{e^x}{e^x-1}\right)dx=\int_0^\infty\frac{x^2e^{-x}}{1-e^{-x}}\ln\left(\frac{1}{1-e^{-x}}\right)dx\\
\overset{e^{-x}=y}{=}-\int_0^1\frac{\ln^2y}{1-y}\ln(1-y)dy=\sum_{n=1}^\infty H_n\int_0^1 y^n \ln^2y dy\\=2\sum_{n=1}^\infty\frac{H_n}{(n+1)^3}=2\left(\sum_{n=1}^\infty\frac{H_n}{n^3}-\sum_{n=1}^\infty \frac1{n^4}\right)\\=2\left(\frac54\zeta(4)-\zeta(4)\right)\\=\frac12\zeta(4)$$
The result $\sum_{n=1}^\infty \frac{H_n}{n^3}=\frac54\zeta(4)$ follows from using the Euler identity:
$$\sum_{n=1}^\infty \frac{H_n}{n^q}= \left(1+\frac{q}{2} \right)\zeta(q+1)-\frac{1}{2}\sum_{k=1}^{q-2}\zeta(k+1)\zeta(q-k)$$

In our calculations, I used the following identities:
$$\sum_{n=1}^\infty H_n x^n=-\frac{\ln(1-x)}{1-x}$$
$$\int_0^1 x^{n-1}\ln^axdx=\frac{(-1)^a a!}{n^{a+1}}$$

Addendum:
The integral $\int_0^1 \frac{\ln^2x}{1-x}\ln(1-x)dx$ can be evaluated without using Euler identity, the first way is by using beta function and the second way is by using the rule
$$\int_0^1\frac{x^{n}\ln^m(x)\ln(1-x)}{1-x}\ dx=\frac12\frac{\partial^m}{\partial n^m}\left(H_n^2+H_n^{(2)}\right)$$
Just set $m=2$ then let $n$ approach $0$ we get
$$\int_0^1 \frac{\ln^2x}{1-x}\ln(1-x)dx=\frac12\frac{\partial^2}{\partial n^2}\left(H_n^2+H_n^{(2)}\right)_{n\to 0}\\=\frac12\left(4H_nH_n^{(3)}+2\left(H_n^{(2)}\right)^2+6H_n^{(4)}-4\zeta(2)H_n^{(2)}-4\zeta(3)H_n-\zeta(4)\right)_{n\to 0}\\=-\frac12\zeta(4)$$
A: More generally,
$\begin{array}\\
\sum_{n=1}^{\infty}\left(\dfrac{1}{n}\sum_{k=n+1}^{\infty}\dfrac{1}{k^m}\right)
&=\sum_{k=2}^{\infty}\sum_{n=1}^{k-1}\dfrac{1}{n}\dfrac{1}{k^m}\\
&=\sum_{k=2}^{\infty}\dfrac{1}{k^m}\sum_{n=1}^{k-1}\dfrac{1}{n}\\
&=\sum_{k=2}^{\infty}\dfrac{1}{k^m}(\sum_{n=1}^{k}\dfrac{1}{n}-\dfrac1{k})\\
&=\sum_{k=2}^{\infty}\dfrac{H_{k}}{k^m}-\sum_{k=2}^{\infty}\dfrac{1}{k^{m+1}}\\
&=\sum_{k=1}^{\infty}\dfrac{H_{k}}{k^m}-\sum_{k=1}^{\infty}\dfrac{1}{k^{m+1}}\\
&=\frac{m+2}{2}\zeta(m+1)-\frac12\sum_{j=1}^{m-2}\zeta(m-j)\zeta(j+1)-\zeta(m+1)
\qquad (*)\\
&=\frac{m}{2}\zeta(m+1)-\frac12\sum_{j=1}^{m-2}\zeta(m-j)\zeta(j+1)\\
\end{array}
$
For $m=3$ this is
$\begin{array}\\
\frac{m}{2}\zeta(m+1)-\frac12\sum_{j=1}^{m-2}\zeta(m-j)\zeta(j+1)
&=\frac{3}{2}\zeta(4)-\frac12\sum_{j=1}^{1}\zeta(3-j)\zeta(j+1)\\
&=\dfrac32\zeta(4)-\frac12\zeta(2)\zeta(2)\\
&=\dfrac32\dfrac{\pi^4}{90}-\frac12(\dfrac{\pi^2}{6})^2\\
&=\pi^4(\dfrac1{60}-\dfrac1{72})\\
&=\pi^4\dfrac{12}{60\cdot 72}\\
&=\pi^4\dfrac{1}{360}\\
\end{array}
$
$(*)$ By the result in this answer:
Generalized Euler sum $\sum_{n=1}^\infty \frac{H_n}{n^q}$
$\sum_{k=1}^\infty\frac{H_k}{k^m}
=\frac{m+2}{2}\zeta(m+1)-\frac12\sum_{j=1}^{m-2}\zeta(m-j)\zeta(j+1)
$
