Infinite Series $\sum\limits_{n=1}^\infty\frac{(H_n)^2}{n^3}$ How to prove that
$$\sum_{n=1}^{\infty}\frac{(H_n)^2}{n^3}=\frac{7}{2}\zeta(5)-\zeta(2)\zeta(3)$$
$H_n$ denotes the harmonic numbers.
 A: One systematic way of proving identities like this is to write everything in terms of multiple harmonic sums and multiple zeta values . For integers $s_1,\ldots,s_k,n\geq 1$, we define the multiple harmonic sum (or MHS)
$$
H_n(s_1,\ldots,s_k):=\sum_{n\geq n_1>\ldots>n_k\geq 1}\frac{1}{n_1^{s_1}\ldots n_k^{s_k}}\in\mathbb{Q}.
$$
The MHS $H_n(1)$ is what you call $H_n$.
We also define the multiple zeta value (or MZV)
$$
\zeta(s_1,\ldots,s_k):=\sum_{n_1>\ldots>n_k\geq 1}\frac{1}{n_1^{s_1}\ldots n_k^{s_k}}\in\mathbb{R},
$$
where we need $s_1\geq 2$ to ensure convergence.
The following relationship between MHS and MZV is easy to check, and will be useful:
$$
\sum_{n=1}^\infty \frac{H_n(s_1,\ldots,s_k)}{n^s}=\zeta(s,s_1,\ldots,s_k)+\zeta(s+s_1,s_2,\ldots,s_k).
$$
MHS's and MZV's satisfy what's called a quasi-shuffle identity. One instance of this is given by the following:
\begin{eqnarray*}
H_n(1)^2&=&\left(\sum_{n_1=1}^n\frac{1}{n_1}\right)\left(\sum_{n_2=1}^n\frac{1}{n_2}\right)\\
&=&\left(\sum_{n\geq n_1>n_2}+\sum_{n\geq n_2>n_1}+\sum_{n\geq n_1=n_2}\right)\frac{1}{n_1n_2}\\
&=&2H_n(1,1)+H_n(2).
\end{eqnarray*}
Using the relationship between MHS and MZV and the quasi-shuffle identity, the left hand side of your identity equals
$$
2\zeta(3,1,1)+2\zeta(4,1)+\zeta(3,2)+\zeta(5).
$$
The expression $\zeta(2)\zeta(3)$ can also be expanded with a quasi-shuffle identity:
$$
\zeta(2)\zeta(3)=\zeta(2,3)+\zeta(3,2)+\zeta(5).
$$
This means your identity is equivalent to the MZV identity
$$
2\zeta(3,1,1)+2\zeta(4,1)+2\zeta(3,2)+\zeta(2,3)=\frac{3}{2}\zeta(5).
$$
This identity is linear and homogeneous (that is, each MZV $\zeta(s_1,\ldots,s_k)$ that appears satisfies $s_1+\ldots+s_k=5$).
There is a ton of literature on producing homogeneous relations among multiple zeta values. One general class of relations, called the extended double shuffle relations, is conjectured to include all relations.
You identity follows by taking a linear combination of the following MZV relations, which can be found in the literature: the double shuffle relation applied to $\zeta(2)\zeta(3)$ (see Derivation and double shuffle relations for multiple zeta values, by Ihara, Kaneko, and Zagier):
$$
\zeta(5)=2\zeta(3,2)+6\zeta(4,1),
$$
duality applied to $\zeta(3,1,1)$ (this is a consequence of an expression of Konstsevich for MZV as iterated integrals, see The algebra of multiple harmonic series, by Hoffman):
$$
\zeta(3,1,1)=\zeta(4,1),
$$
and the sum formula (known in this case by Euler and proven in the general case independently by Granville and Zagier, the statement can also be found in The algebra of multiple harmonic series):
$$
\zeta(2,3)+\zeta(3,2)+\zeta(4,1)=\zeta(5).
$$
A: This is a known result; an explicit reference is [PP, equation (3c)].
The proof there is elementary (no need for digamma functions, 
definite integrals and contour integrals, etc.), but nontrivial, requiring 
clever manipulations with identities such as $1/XY = 1/(X(X+Y)) + 1/(Y(X+Y))$.
As Julian Rosen noted, such sums are often expressed in terms of
multiple
zeta functions, such as the double and triple zetas
$$
\zeta(a,b) = \mathop{\sum\sum}_{0<m<n} \frac1{m^a n^b},
\quad
\zeta(a,b,c) = \mathop{\sum\sum\sum}_{0<l<m<n} \frac1{l^a m^b n^c}.
$$
I found [PP] via Michael Hoffmann's list of

References on multiple zeta values and Euler sums.
[PP] cites a paper of Borwein and Girgensohn [BG],
where the key triple-zeta value
$\zeta(3,1,1) = 2 \zeta(5) - \zeta(3) \zeta(2)$
is given on page 21, together with the note that
all such values of weight at most 6 appear in [M].
Indeed here the weight is $3+1+1 = 5 \leq 6$ and the result is
the case $p=3$ of Theorem 4.1, listed explicitly as
equation (4.2) on page 126.
References
[BG]
J. Borwein and R. Girgensohn,

Evaluation of triple Euler sums,
Electronic Journal of Combinatorics 3, research paper #23, 1996.
[M] C. Markett:

Triple sums and the Riemann zeta function,
J. Number Theory 48 (1994), 113$-$132.
[PP] Alois Panholzer and Helmut Prodinger:

Computer-free evaluation of an infinite double sum via Euler sums,
Séminaire Lotharingien de Combinatoire 55 (2005), Article B55a.
A: using the following identity :
$$\displaystyle \frac{\ln^2(1-x)}{1-x}=\sum_{n=1}^{\infty}\left(H_n^2-H_n^{(2)}\right)x^n$$
multiply both sides by $\frac{\ln^2x}{x}$ then integrate both sides w.r.t $x$ from $0$ to $1$, we have:
\begin{align*}
S&=\sum_{n=1}^{\infty}\left(H_n^2-H_n^{(2)}\right)\int_0^1x^{n-1}\ln^2x\ dx=\color{blue}{2\sum_{n=1}^{\infty}\frac{H_n^2-H_n^{(2)}}{n^3}}=\int_0^1\frac{\ln^2x\ln^2(1-x)}{x(1-x)}\ dx\\
&=\int_0^1\frac{\ln^2x\ln^2(1-x)}{x}\ dx+\underbrace{\int_0^1\frac{\ln^2x\ln^2(1-x)}{1-x}\ dx}_{x\mapsto1-x}=2\int_0^1\frac{\ln^2x\ln^2(1-x)}{x}\ dx\\
&=4\sum_{n=1}^{\infty}\left(\frac{H_n}{n}-\frac1{n^2}\right)\int_0^1x^{n-1}\ln^2x\ dx=4\sum_{n=1}^{\infty}\left(\frac{H_n}{n}-\frac1{n^2}\right)\frac{2}{n^3}=\color{blue}{8\sum_{n=1}^{\infty}\frac{H_n}{n^4}-8\zeta(5)}
\end{align*}
where we used $\ln^2(1-x)=2\sum_{n=1}^\infty\frac{H_n}{n+1}x^{n+1}=2\sum_{n=1}^\infty\left(\frac{H_n}{n}-\frac1{n^2}\right)x^n$
rearranging the terms of the blue sides, we have:
\begin{align*}
\sum_{n=1}^{\infty}\frac{H_n^2}{n^3}&=\sum_{n=1}^{\infty}\frac{H_n^{(2)}}{n^3}+4\sum_{n=1}^{\infty}\frac{H_n}{n^4}-4\zeta(5)\\
&=\left(3\zeta(2)\zeta(3)-\frac92\zeta(5)\right)+4\left(3\zeta(5)-\zeta(2)\zeta(3)\right)-4\zeta(5)\\
&=\frac72\zeta(5)-\zeta(2)\zeta(3)
\end{align*}
The proof of $\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^3}=3\zeta(2)\zeta(3)-\frac92\zeta(5)$ can be found here.
A: Since $$\sum_{n=1}^\infty\frac{H_{n-1}}{n}x^n=\frac12\ln^2(1-x)$$
divide both sides by $x$ then $\int_0^y$
$$\sum_{n=1}^\infty\frac{H_{n-1}}{n^2}y^n=\frac12\int_0^y\frac{\ln^2(1-x)}{x}dx$$
Next multiply both sides by $-\frac{\ln(1-y)}{y}$ then $\int_0^1$ and apply that $-\int_0^1 y^{n-1}\ln(1-y)dy=\frac{H_n}{n}$
$$\sum_{n=1}^\infty\frac{H_{n-1}H_n}{n^3}=\sum_{n=1}^\infty\frac{H_n^2}{n^3}-\sum_{n=1}^\infty\frac{H_n}{n^4}=-\frac12\int_0^1\int_0^y\frac{\ln^2(1-x)\ln(1-y)}{xy}dxdy$$
$$=-\frac12\int_0^1\frac{\ln^2(1-x)}{x}\left(\int_x^1\frac{\ln(1-y)}{y}dy\right)dx$$
$$=-\frac12\int_0^1\frac{\ln^2(1-x)}{x}\left(\text{Li}_2(x)-\zeta(2)\right)dx$$
$$=\zeta(2)\zeta(3)-\frac12\int_0^1\frac{\ln^2(1-x)\text{Li}_2(x)}{x}dx$$
$$=\zeta(2)\zeta(3)-\frac12\sum_{n=1}^\infty\frac1{n^2}\int_0^1 x^{n-1}\ln^2(1-x)dx$$
$$=\zeta(2)\zeta(3)-\frac12\sum_{n=1}^\infty\frac1{n^2}\left(\frac{H_n^2+H_n^{(2)}}{n}\right)$$
$$=\zeta(2)\zeta(3)-\frac12\sum_{n=1}^\infty\frac{H_n^2}{n^3}-\frac12\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^3}$$
Rearrange the terms
$$\sum_{n=1}^\infty\frac{H_n^2}{n^3}=\frac23\zeta(2)\zeta(3)+\frac23\sum_{n=1}^\infty\frac{H_n}{n^4}-\frac13\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^3}=\frac{7}{2}\zeta(5)-\zeta(2)\zeta(3)$$
where we used $\sum_{n=1}^\infty\frac{H_n}{n^4}=3\zeta(5)-\zeta(2)\zeta(3)$ and $\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^3}=3\zeta(2)\zeta(3)-\frac92\zeta(5)$
