Infinite Series $\sum\limits_{n=1}^\infty\left(\frac{H_n}n\right)^2$ How can I find a closed form for the following sum?
$$\sum_{n=1}^{\infty}\left(\frac{H_n}{n}\right)^2$$
($H_n=\sum_{k=1}^n\frac{1}{k}$).
 A: Starting with $\displaystyle \dfrac{H_n}{n} = \sum_{k=1}^{\infty} \dfrac{1}{k(k+n)}$ we have,
\begin{align*}\sum_{n=1}^{\infty} \dfrac{H_n^2}{n^2} &= \sum_{n=1}^{\infty} \left(\sum_{k=1}^{\infty}\dfrac{1}{k(k+n)}\right)^2\\&= \sum_{n=1}^{\infty}\sum_{k,j=1}^{\infty} \dfrac{1}{jk(j+n)(k+n)} \\&= \sum_{n=1}^{\infty} \left(\sum_{k=1}^{\infty} \dfrac{1}{k^2(n+k)^2}+ 2\sum_{1 \le k<j} \dfrac{1}{jk(j+n)(k+n)}\right) \\&= \sum_{1 \le k < j} \dfrac{1}{k^2j^2} + 2\sum_{n=1}^{\infty}\sum_{k,m=1}^{\infty} \dfrac{1}{k(k+m)(k+m+n)(k+n)}\\&= \dfrac{1}{2}\left(\left(\sum_{k=1}^{\infty}\dfrac{1}{k^2}\right)^2 – \sum_{k=1}^{\infty}\dfrac{1}{k^4}\right) + 2\sum_{k,m,n=1}^{\infty}\dfrac{(k+m)(k+n) – k(k+m+n)}{kmn(k+m)(k+m+n)(k+n)} \\&= \dfrac{1}{2}\left(\zeta^2(2) – \zeta(4)\right) + 2\sum_{k,m,n=1}^{\infty} \dfrac{1}{kmn(k+m+n)} – 2\sum_{k,m,n=1}^{\infty} \dfrac{1}{mn(m+k)(n+k)}\end{align*}
Therefore, $\displaystyle 3\sum_{n=1}^{\infty} \left(\sum_{k=1}^{\infty}\dfrac{1}{k(k+n)}\right)^2 = \dfrac{1}{2}\left(\zeta^2(2) - \zeta(4)\right) + 2\sum_{k,m,n=1}^{\infty} \dfrac{1}{kmn(k+m+n)}$
Using, $\displaystyle \sum_{k,m,n=1}^{\infty} \dfrac{1}{kmn(k+m+n)} = 6\zeta(4)$ we conclude, 

$$\sum_{n=1}^{\infty} \dfrac{H_n^2}{n^2} = \dfrac{1}{6}\zeta^2(2) + \dfrac{23}{6}\zeta(4)$$

To see the last result, \begin{align*} \sum_{k,m,n=1}^{\infty} \dfrac{1}{kmn(k+m+n)} &= \sum_{k,m,n=1}^{\infty} \int_0^1 \dfrac{x^{k+m+n}}{kmn}\,\dfrac{dx}{x} \\&= \int_0^1 \log^3(1-x)\,\dfrac{dx}{x} \\&= -\int_0^1 \dfrac{\log^3 x}{1-x}\,dx \\&= -\sum_{n=0}^{\infty} \int_0^1 x^n\log^3 x\,dx \\&= 6\sum_{n=0}^{\infty}\dfrac{1}{(n+1)^4} = 6\zeta(4)\end{align*}
A: I think this is the shortest solution 
Using the generating function
$$\frac12\ln^2(1-x)=\sum_{n=1}^\infty\frac{H_n}{n+1}x^{n+1}=\sum_{n=1}^\infty\frac{H_{n-1}}{n}x^n$$
Multiply both sides by $\frac{\ln(1-x)}{x}$ then $\int_0^1$ and use the fact that $\int_0^1 x^{n-1}\ln(1-x)=-\frac{H_n}{n}$ we get
$$\frac12\int_0^1\frac{\ln^3(1-x)}{x}\ dx=-3\zeta(4)=\sum_{n=1}^\infty\frac{H_{n-1}}{n}\left(-\frac{H_n}{n}\right)=\sum_{n=1}^\infty\frac{H_n}{n^3}-\sum_{n=1}^\infty\frac{H_n^2}{n^2}$$
Substituting $\displaystyle\sum_{n=1}^\infty\frac{H_n}{n^3}=\frac54\zeta(4)$ gives $\displaystyle\sum_{n=1}^\infty\frac{H_n^2}{n^2}=\frac{17}4\zeta(4)$
Note that the fuction used above follows from integrating both sides of $\sum_{n=1}^\infty x^n H_n=-\frac{\ln(1-x)}{1-x}$
A: EDITED. Some simplifications were made.

Here is a solution.
1. Basic facts on the dilogarithm. Let $\mathrm{Li}_{2}(z)$ be the dilogarithm function defined by
$$ \operatorname{Li}_{2}(z) = \sum_{n=1}^{\infty} \frac{z^{n}}{n^{2}} = - \int_{0}^{z} \frac{\log(1-x)}{x} \, dx. $$
Here the branch cut of $\log $ is chosen to be $(-\infty, 0]$ so that $\operatorname{Li}_{2}$ defines a holomorphic function on the region $\Bbb{C} \setminus [1, \infty)$. Also, it is easy to check (by differentiating both sides) that the following identities hold
\begin{align*}
\operatorname{Li}_{2}\left(\tfrac{z}{z-1}\right)
&= -\mathrm{Li}_{2}(z) - \tfrac{1}{2}\log^{2}(1-z); \quad z \notin [1, \infty) \tag{1} \\
\operatorname{Li}_{2}\left(\tfrac{1}{1-z}\right)
&= \color{blue}{\boxed{\operatorname{Li}_{2}(z) + \zeta(2) - \tfrac{1}{2}\log^{2}(1-z)}} + \color{red}{\boxed{\log(-z)\log(1-z)}}; \quad z \notin [0, \infty) \tag{2}
\end{align*}
Notice that in (2), the blue-colored part is holomorphic on $|z| < 1$ while the red-colored part induces the branch cut $[-1, 0]$.
2. A useful power series. Now let us consider the power series
$$ f(z) = \sum_{n=0}^{\infty} \frac{H_n}{n} z^n. $$
Then $f(z)$ is automatically holomorphic inside the disc $|z| < 1$. Moreover, it is easy to check that
$$ \sum_{n=1}^{\infty} H_{n} z^{n-1}
= \frac{1}{z} \left( \sum_{n=1}^{\infty} \frac{z^{n}}{n} \right)\left( \sum_{n=0}^{\infty} z^{n}\right)
= -\frac{\log(1-z)}{z(1-z)}. $$
thus integrating both sides, together with the identity $\text{(1)}$, we obtain the following representation of $f(z)$.
$$f(z)
= \operatorname{Li}_{2}(z) + \tfrac{1}{2}\log^{2}(1-z)
= -\operatorname{Li}_{2}\left(\tfrac{z}{z-1}\right). \tag{3}$$
3. Integral representation and the result. By the Parseval's identity, we have
$$ \sum_{n=1}^{\infty} \frac{H_{n}^{2}}{n^{2}}
= \frac{1}{2\pi} \int_{0}^{2\pi} f(e^{it})f(e^{-it}) \, dt
= \frac{1}{2\pi i} \int_{|z|=1} \frac{f(z)}{z} f\left(\frac{1}{z}\right) \, dz \tag{4} $$
Since $\frac{1}{z}f(z)$ is holomorphic inside $|z| = 1$, the failure of holomorphy of the integrand stems from the branch cut of
\begin{align*}
f\left(\tfrac{1}{z}\right)
&= -\operatorname{Li}_{2}\left(\tfrac{1}{1-z}\right) \\
&= -\color{blue}{\left( \operatorname{Li}_{2}(z) + \zeta(2) - \tfrac{1}{2}\log^{2}(1-z) \right)} - \color{red}{\log(-z)\log(1-z)},
\end{align*}
which is $[0, 1]$. To resolve this, we utilize the identity $\text{(2)}$. Note that the blue-colored portion does not contributes to the the integral $\text{(4)}$, since it remains holomorphic inside $|z| < 1$. That is, only the red-colored portion gives contribution to the integral. Consequently we have
\begin{align*}
\sum_{n=1}^{\infty} \frac{H_{n}^{2}}{n^{2}}
&= -\frac{1}{2\pi i} \int_{|z|=1} \frac{f(z)}{z} \color{red}{\log(-z)\log(1-z)} \, dz. \tag{5}
\end{align*}
Since the integrand is holomorphic on $\Bbb{C} \setminus [0, \infty)$, we can utilize the keyhole contour wrapping around $[0, 1]$ to reduce $\text{(5)}$ to
\begin{align*}
\sum_{n=1}^{\infty} \frac{H_{n}^{2}}{n^{2}}
&=-\frac{1}{2\pi i} \Bigg\{ \int_{0^{-}i}^{1+0^{-}i} \frac{f(z)\log(-z)\log(1-z)}{z} \, dz \\
&\qquad \qquad + \int_{1+0^{+}i}^{+0^{+}i} \frac{f(z)\log(-z)\log(1-z)}{z} \, dz \Bigg\} \\
&=-\frac{1}{2\pi i} \Bigg\{ \int_{0}^{1} \frac{f(x)(\log x + i\pi)\log(1-x)}{x} \, dx \\
&\qquad \qquad - \int_{0}^{1} \frac{f(x)(\log x - i\pi)\log(1-x)}{x} \, dx \Bigg\} \\
&=-\int_{0}^{1} \frac{f(x)\log(1-x)}{x} \, dx. \tag{5}
\end{align*}
Plugging $\text{(3)}$ to the last integral and simplifying a little bit, we have
\begin{align*}
\sum_{n=1}^{\infty} \frac{H_{n}^{2}}{n^{2}}
&= - \int_{0}^{1} \frac{\operatorname{Li}_2(x)\log(1-x)}{x} \, dx - \frac{1}{2}\int_{0}^{1} \frac{\log^{3}(1-x)}{x} \, dx \\
&= \left[ \frac{1}{2}\operatorname{Li}_2(x)^2 \right]_0^1 - \frac{1}{2} \int_{0}^{1} \frac{\log^3 x}{1-x} \, dx \\
&= \frac{1}{2}\zeta(2)^{2} + \frac{1}{2} \Gamma(4)\zeta(4) \\
&= \frac{17\pi^{4}}{360}
\end{align*}
as desired.
A: Different approach:
Start with the identity
$$\sum_{n=1}^\infty (H_n^{(2)}-H_n^2)x^{n}=-\frac{\ln^2(1-x)}{1-x}$$
Multiply both sides by $-\frac{\ln x}{x}$ and integrate between $0$ and $1$ and use $\int_0^1-x^{n-1}\ln x\ dx=\frac1{n^2}$ we get
$$\sum_{n=1}^\infty \frac{H_n^{(2)}-H_n^2}{n^2}=\int_0^1\frac{\ln x\ln^2(1-x)}{x(1-x)}dx=\int_0^1\frac{\ln(1-x)\ln^2x}{(1-x)x}dx$$
$$=-\sum_{n=1}^\infty H_n\int_0^1 x^{n-1}\ln^2x\ dx=-2\sum_{n=1}^\infty\frac{H_n}{n^3}=-\frac52\zeta(4)$$
$$\Longrightarrow\sum_{n=1}^\infty\frac{H_n^2}{n^2}=\frac52\zeta(4)+\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^2}=\frac52\zeta(2)+\frac74\zeta(4)=\frac{17}4\zeta(4)$$
A: This result was a problem in American mathematical monthly in 50s.
Source:
H. F. Sandham and Martin Kneser, The American mathematical monthly, Advanced problem 4305, Vol. 57, No. 4 (Apr., 1950), pp. 267-268
\begin{align}S&=\sum_{n=1}^\infty \left(\frac{\text{H}_n}{n}\right)^2\\
&=\sum_{n=1}^\infty \frac{1}{n^2}\left(\int_0^1 \frac{1-t^n}{1-t}dt\right)\left(\int_0^1 \frac{1-u^n}{1-u}du\right)\\
&=\int_0^1 \int_0^1 \frac{\text{Li}_2(1)+\text{Li}_2( tu)-\text{Li}_2(t)-\text{Li}_2(u)}{(1-t)(1-u)}dtdu\\
&\overset{\text{IBP}}=\int_0^1 \int_0^1 \frac{\ln(1-t)\big(\ln(1-t)-\ln(1-tu)\big)}{t(1-u)}dtdu\\ 
&\overset{x=1-tu,y=\frac{1-t}{1-tu}}=\int_0^1\int_0^1 \frac{\ln y\ln(xy)}{(1-y)(1-xy)}dxdy\\
&\overset{w(x)=yx}=\int_0^1 \frac{\ln y}{y(1-y)}\left(\int_0^y \frac{\ln w}{1-w}dw\right)dy\\
&=\int_0^1 \frac{\ln y}{1-y}\left(\int_0^y \frac{\ln w}{1-w}dw\right)dy+\underbrace{\int_0^1 \frac{\ln y}{y}\left(\int_0^y \frac{\ln w}{1-w}dw\right)dy}_{\text{IBP}}\\
&=\frac{1}{2}\left(\int_0^1 \frac{\ln w}{1-w}dw\right)^2-\frac{1}{2}\int_0^1 \frac{\ln^3 w}{1-w}dw\\
&=\frac{1}{2}\times \left(\frac{\pi^2}{6}\right)^2-\frac{1}{2}\times -\frac{\pi^4}{15}\\
&=\boxed{\dfrac{17\pi^4}{360}}.
\end{align} NB:
I assume that:
\begin{align}\zeta(2)=\frac{\pi^2}{6},&&\zeta(4)=\frac{\pi^4}{90},&&\int_0^1 \frac{\ln w}{1-w}dw=-\zeta(2)=-\frac{\pi^2}{6},&&\int_0^1 \frac{\ln^3 w}{1-w}dw=-6\zeta(4)=-\frac{\pi^4}{15}\end{align}
A: From here we have  
$$\displaystyle\int_0^1 x^{n-1}\ln^2(1-x)\ dx=\frac1n\left({H_n^2}+H_n^{(2)}\right)$$
dividing both sides by $n$ then summing w.r.t $n$ from $n=1$ to $\infty$ we get
\begin{align*}
\sum_{n=1}^{\infty}\frac1{n^2}\left({H_n^2}+H_n^{(2)}\right)&=\int_0^1\frac{\ln^2(1-x)}{x}\sum_{n=1}^{\infty}\frac{x^n}{n}\ dx=-\int_0^1\frac{\ln^3(1-x)}{x}\ dx\\
&=-\int_0^1\frac{\ln^3(x)}{1-x}\ dx=6\sum_{n=1}^{\infty}\frac{1}{n^4}=6\zeta(4)
\end{align*}
we have, using $\displaystyle\sum_{n=1}^{\infty}\frac{H_n^{(a)}}{n^a}=\frac12\left(\zeta(2a)+\zeta^2(a)\right)$ that  $\displaystyle\sum_{n=1}^{\infty}\frac{H_n^{(2)}}{n^2}=\frac12\left(\zeta(4)+\zeta^2(2)\right)=\frac74\zeta(4)$
finally $$\displaystyle\sum_{n=1}^{\infty}\frac{H_n^2}{n^2}=6\zeta(4)-\frac74\zeta(4)=\frac{17}4\zeta(4)$$
A: I believe the answer you're looking for is in this Wikipedia article :

The following identity was first conjectured by Enrico Au-Yeung , a student of Jonathan Borwein, using computer search and the PSLQ algorithm, in 1993 :
  $$\sum_{k=1}^\infty \frac{1}{k^2}\left(1+\frac{1}{2}+\frac{1}{3}+\cdots+\frac{1}{k}\right)^2 = \frac{17\pi^4}{360}.$$

A simple Google search will return several papers in PDF format containing this and other curious and interesting mathematical identities. Or you could simply visit David H. Bailey's own page, and search for papers containing the string experiment in their title, most of which also contain this and many other similar results. The proofs are based on a combination of one or more of the following: the PSLQ algorithm I've already mentioned, computer-assisted proofs, and-or inverse symbolic computation.
A: Using the generating function:
\begin{gather}
\sum_{n=1}^\infty\frac{H_{n}^2}{n^2}x^{n}=\operatorname{Li}_4(x)-2\operatorname{Li}_4(1-x)+2\ln(1-x)\operatorname{Li}_3(1-x)+\frac12\operatorname{Li}_2^2(x)\nonumber\\
\qquad-\ln^2(1-x)\operatorname{Li}_2(1-x)-\frac13\ln(x)\ln^3(1-x)+2\zeta(4),
\end{gather}
we have the following results:
$$\sum_{n=1}^\infty \frac{H_n^{2}}{n^2}=\frac{17}{4}\zeta(4);$$
$$\sum_{n=1}^{\infty}\frac{(-1)^nH_n^2}{n^2}=2\operatorname{Li}_4\left(\frac12\right)-\frac{41}{16}\zeta(4)+\frac74\ln(2)\zeta(3)-\frac12\ln^2(2)\zeta(2)+\frac1{12}\ln^4(2);$$
$$\sum_{n=1}^\infty\frac{H_n^2}{n^22^n}=-\operatorname{Li}_4\left(\frac12\right)+\frac{37}{16}\zeta(4)-\frac74\ln(2)\zeta(3)+\frac14\ln^2(2)\zeta(2)-\frac1{24}\ln^4(2).$$
