Let $H_{n}$ be the nth harmonic number defined by $ H_{n} := \sum_{k=1}^{n} \frac{1}{k}$.
How would you prove that
$$\sum_{n=1}^{\infty}\frac{H_n}{n^3}=\frac{\pi^4}{72}?$$
Simply replacing $H_{n}$ with $\sum_{k=1}^{n} \frac{1}{k}$ does not seem like a good starting point. Perhaps another representation of the nth harmonic number would be more useful.
The Euler sum $\sum_{n=1}^{\infty} \frac{H_{n}}{n^{q}}$, where $q$ is an odd positive integer greater than $1$, can also be evaluated using this approach. See here.
Using the integration representation $$H_{n} = \int_{0}^{1} \frac{1-t^{n}}{1-t} \, dt \ ,$$
we have
$$ \begin{align} \sum_{n=1}^{\infty} \frac{H_{n}}{n^{3}} &= \sum_{n=1}^{\infty} \frac{1}{n^{3}} \int_{0}^{1} \frac{1-t^{n}}{1-t} \ dt \\ &= \int_{0}^{1} \frac{1}{1-t} \sum_{n=1}^{\infty} \frac{1-t^{n}}{n^{3}} \\ &=\int_{0}^{1} \frac{\zeta(3) - \text{Li}_{3}(t)}{1-t} \ dt \tag{1} \\ &=- \Big(\zeta(3)-\text{Li}_{3}(t)\Big) \ln(1-t) \Bigg|_{0}^{1} - \int^{1}_{0} \frac{ \text{Li}_{2}(t) \log(1-t)}{t} \ dt \\ &= -\int_{0}^{1} \frac{ \text{Li}_{2}(t) \log(1-t)}{t} \ dt \\ &= \int_{0}^{1} \text{Li}_{2}(t) \, d \big(\text{Li}_{2}(t)\big) \\ &= \frac{\big(\text{Li}_{2}(t)\big)^{2}}{2} \Bigg|^{1}_{0} \\ &= \frac{\zeta^{2}(2)}{2} \\ &= \frac{\pi^{4}}{72}. \end{align}$$
$ $
I will try to reduce the sum to an integral: $$ \sum_{n=1}^\infty \frac{H_n}{n^3} = \sum_{n=1}^\infty H_n \frac{1}{\Gamma(3)} \int_0^\infty x^2 \mathrm{e}^{-n x} \mathrm{d} x = \frac{1}{2} \int_0^\infty x^2 \sum_{n=1}^\infty H_n \mathrm{e}^{-n x} \mathrm{d} x \tag{1} $$ We now make use of the following generating function: $$ \sum_{n=1}^\infty H_n z^n = \sum_{n=1}^\infty H_n \Delta_n \left(\frac{z^n}{z-1} \right) $$ where $\Delta_n f_n = f_{n+1}-f_n$. We can now use summation by parts: $$ \sum_{n=1}^m a_n \Delta_n b_n = b_{m+1} a_m - b_1 a_1 - \sum_{n=1}^{m-1} b_{n+1} \Delta_n a_n $$ with $b_n = \frac{z^n}{z-1}$ and $a_n = H_n$, and using $\Delta_n H_n = \frac{1}{n+1}$, we get $$ \sum_{n=1}^\infty H_n z^n = \sum_{n=1}^\infty H_n \Delta_n \left(\frac{z^n}{z-1} \right) = -1 - \sum_{n=1}^\infty \frac{z^{n+1}}{z-1} \frac{1}{n+1} = \frac{\log(1-z)}{z-1} \tag{2} $$ Now, using $(2)$ in $(1)$: $$ \sum_{n=1}^\infty \frac{H_n}{n^3} = -\frac{1}{2} \int_0^\infty x^2 \frac{\log\left(1-\mathrm{e}^{-x}\right)}{1-\mathrm{e}^{-x}} \mathrm{d}x \stackrel{t=\exp(-x)}{=} -\frac{1}{2} \int_0^1 \frac{\log(1-t)}{1-t} \frac{\log^2(t)}{t} \mathrm{d}t \tag{3} $$ The latter integral can be evaluated using derivatives of the Euler beta function: $$ \int_0^1 \frac{\log(1-t)}{1-t} \frac{\log^2(t)}{t} \mathrm{d}t = \lim_{\alpha \downarrow 0} \lim_{\beta \downarrow 0} \frac{\mathrm{d}}{\mathrm{d} \alpha} \frac{\mathrm{d}^2}{\mathrm{d} \beta^2} \int_0^1 \left(1-t\right)^{\alpha-1} t^{\beta-1} \mathrm{d} t = \lim_{\alpha \downarrow 0} \lim_{\beta \downarrow 0} \frac{\mathrm{d}}{\mathrm{d} \alpha} \frac{\mathrm{d}^2}{\mathrm{d} \beta^2} \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)} $$ Using $$ \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)} = \left(\frac{1}{\alpha} + \frac{1}{\beta} \right) \frac{\Gamma(\alpha+1) \Gamma(\beta+1)}{\Gamma(\alpha + \beta+1)} = \left(\frac{1}{\alpha} + \frac{1}{\beta} \right) \left( 1 - \frac{\pi^2}{6} \alpha \beta + \left(\alpha \beta^2 + \beta \alpha^2\right) \zeta(3) - \frac{\pi^4}{360} \left(4 \alpha \beta^3 + \alpha^2 \beta^2 + 4 \alpha^3 \beta\right) + \cdots \right) $$ Differentiating we get the result: $$ \lim_{\alpha \downarrow 0} \lim_{\beta \downarrow 0} \frac{\mathrm{d}}{\mathrm{d} \alpha} \frac{\mathrm{d}^2}{\mathrm{d} \beta^2} \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)} = -\frac{\pi^2}{36} $$ yielding with eq. $(3)$: $$ \sum_{n=1}^\infty \frac{H_n}{n^3} = \frac{\pi^4}{72} $$
Apparently Euler showed in 1775 that: $$2 \sum_{n=1}^{\infty}\frac{H_n}{n^q} = (q+2)\zeta(q+1)- \sum_{m=1}^{q-2}\zeta(m+1)\zeta(q-m)$$ In your case, $q=3$, so that: $$2 \sum_{n=1}^{\infty}\frac{H_n}{n^3} = 5\zeta(4)- \zeta(2)^2 = 5\frac{\pi^4}{90}-\frac{\pi^4}{36}=\frac{\pi^4}{36}$$ The original proof of euler's in english can be found here.
$\newcommand{\+}{^{\dagger}} \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert} \newcommand{\wt}[1]{\widetilde{#1}}$ $\ds{\sum_{n = 1}^{\infty}{H_{n} \over n^{3}} = {\pi^{4} \over 72}:\ {\large ?}}$
\begin{align}\color{#66f}{\large\sum_{n = 1}^{\infty}{H_{n} \over n^{3}}} &=\sum_{n = 1}^{\infty}{1 \over n^{3}}\ \overbrace{\quad\sum_{k = 1}^{\infty}\pars{{1 \over k} - {1 \over k + n}}\quad} ^{\ds{=\ H_{n}}}\ =\ \sum_{n = 1}^{\infty}\sum_{k = 1}^{\infty}{1 \over n^{2}k\pars{k + n}} \\[3mm]&=\half\sum_{n = 1}^{\infty}\sum_{k = 1}^{\infty} \bracks{{1 \over n^{2}k\pars{k + n}} + {1 \over k^{2}n\pars{n + k}}} =\half\sum_{n = 1}^{\infty}\sum_{k = 1}^{\infty} {k + n \over n^{2}k^{2}\pars{k + n}} \\[3mm]&=\half\sum_{n = 1}^{\infty}\sum_{k = 1}^{\infty}{1 \over n^{2}k^{2}} =\half\pars{\sum_{n = 1}^{\infty}{1 \over n^{2}}}^{2} =\half\pars{\pi^{2} \over 6}^{2}=\color{#66f}{\Large{\pi^{4} \over 72}} \approx 1.3529 \end{align}
Appealing to the integral identity I introduced in a previous problem for the more general case of your sum
$$ B(p,q) = \sum_{k=1}^{\infty} \dfrac{H_k^{(p)}}{k^q}=\frac{(-1)^q}{\Gamma(q)}\int_{0}^{1}\!{\frac {\left(\ln\left(u\right)\right)^{q-1}{Li_{p}(u)} }{ u\left( u-1 \right)}}{du}, $$
where $Li_{p}(u)$ is the polylogarithm function, we can have the following integral representation for the sum
$$B(1,3) = \sum_{k=1}^{\infty} \dfrac{H_k}{k^3}=-\frac{1}{\Gamma(3)}\int_{0}^{1}\!{\frac {\ln^2\left(u\right){\ln(1-u)} }{ u\left( 1-u \right)}}{du} $$
$$ = -\frac{1}{\Gamma(3)}\lim_{w\to 0}\lim_{s \to 0}\frac{d}{dw} \frac{d^2}{ds^2}\int_{0}^{1} u^{s-1}\, (1-u)^{w-1} $$
$$ = -\frac{1}{2}\lim_{w\to 0}\lim_{s \to 0}\frac{d}{dw} \frac{d^2}{ds^2}\beta(s, w)=\frac{\pi^4}{72}, $$
where $\beta(s,w)=\frac{\Gamma(s)\Gamma(w)}{\Gamma(s+w)}$ is the beta function.
Note: $$ Li_{1}(x) = -\ln(1-x). $$
We proved here
$$\sum_{k=1}^\infty\frac{H_k}{k^n}=\frac12\sum_{i=1}^{n-2}(-1)^{i-1}\zeta(n-i)\zeta(i+1),\quad n=3,5,7, ...$$
Set $n=3$ we get $$\sum_{k=1}^\infty \frac{H_k}{k^3}=\frac54\zeta(4)$$
where $\zeta(4)=\frac{\pi^4}{90}$
