# 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}$).

• why would there be a closed form? – Will Jagy Nov 6 '13 at 7:02
• Mathematica says $\frac{17\pi^4}{360}$...good luck figuring that one out. – Tim Ratigan Nov 6 '13 at 7:03
• @Tim.Ratigan Depressing but honest/truthful comment. – Lord Soth Nov 6 '13 at 7:24
• – qoqosz Nov 6 '13 at 8:10
• (harmonic-numbers) seems to be quite reasonable tag for this question. However, there is a limit of at most 5 tags. – Martin Sleziak Aug 21 '14 at 7:05

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.

• $\beta\left(0,1\right)$ diverges. – Felix Marin Jul 27 '14 at 5:29

SOS always has the most clever and ingenious solutions, but if I may contribute something I have found interesting. A fun method of evaluating a whole slew of Euler sums is by using the residues of digamma.

By noting the identity, $\displaystyle \sum_{n=1}^{\infty}\frac{(H_{n})^{2}}{n^{2}}=2\sum_{n=1}^{\infty}\frac{H_{n}}{n^{3}}+\sum_{n=1}^{\infty}\frac{H_{n}^{(2)}}{n^{2}}......[1]$

one can evaluate each of the sums on the right side and thus arrive at the quadratic Euler sum in question.

For the first sum on the right, begin by considering $\displaystyle f(z)=\frac{\left(\gamma+\psi(-z)\right)^{2}}{z^{2}}$ and, due to the poles of digamma, compute the residue at n (the positive integers).

As $z\to n$, the series is $\displaystyle\frac{1}{(z-n)^{2}}+\frac{2H_{n}}{z-n}+\cdot\cdot\cdot$

Thus, the residue is $\displaystyle\lim_{z\to n}\left[Res\left(\frac{1}{(z-n)^{2}}\cdot \frac{1}{z^{3}}\right)+Res\left(\frac{2H_{n}}{z-n}\cdot \frac{1}{z^{3}}\right)\right]$

$\displaystyle=\frac{-3}{n^{4}}+\frac{2H_{n}}{n^{3}}$

Sum these residues: $\displaystyle-3\sum_{n=1}^{\infty}\frac{1}{n^{4}}+2\sum_{n=1}^{\infty}\frac{2H_{n}}{n^{2}}$

By taking the Laurent expansion of f(z), the residue at z=0 is the coefficient of the 1/z term.

$\displaystyle \psi(-z)+\gamma = \frac{1}{z}-\zeta(2)z+\zeta(3)z^{2}-\zeta(4)z^{3}+\cdot\cdot\cdot$

$\displaystyle f(z)=\frac{1}{z^{5}}-\frac{\pi^{2}}{3}\cdot \frac{1}{z^{3}}-2\zeta(3)\cdot \frac{1}{z^{2}}+\frac{\pi^{4}}{180}\cdot \frac{1}{z}+\cdot\cdot\cdot$

As can be seen, the residue at 0 is $\frac{\pi^{4}}{180}$

Put them together, set to 0, and get

$\displaystyle2\sum_{n=1}^{\infty}\frac{2H_{n}}{n^{3}}-3\sum_{n=1}^{\infty}\frac{1}{n^{4}}+\frac{\pi^{4}}{180}=0$

$\displaystyle2H-\frac{\pi^{4}}{30}+\frac{\pi^{4}}{180}=0$

$\displaystyle \boxed{\displaystyle\sum_{n=1}^{\infty}\frac{H_{n}}{n^{3}}=\frac{\pi^{4}}{72}}.......[2]$

Now, for the other sum on the right of [1], where $\displaystyle H_{n}^{(2)}=\sum_{k=1}^{n}\frac{1}{k^{2}}$

$\displaystyle \sum_{n=1}^{\infty}\frac{H_{n}^{(2)}}{n^{2}}$

Due to symmetry of Euler sums, if we have a sum $\displaystyle S_{p,q}=\sum_{n=1}^{\infty}\frac{H_{n}^{(p)}}{n^{q}}$, and $p=q$, then by symmetry $S_{p,q}+S_{q,p}=\zeta(p)\zeta(q)+\zeta(p+q)$

So, in this case with $p=q=2$, then

$\displaystyle2\sum_{n=1}^{\infty}\frac{H_{n}^{(2)}}{n^{2}}=\frac{\pi^{2}}{36}+\frac{\pi^{4}}{90}=\frac{7\pi^{4}}{180}$

$\displaystyle \boxed{\displaystyle\sum_{n=1}^{\infty}\frac{H_{n}^{(2)}}{n^{2}}=\frac{7\pi^{4}}{360}}$

Now, add this to the result of the other sum in [2]:

$\displaystyle\frac{7\pi^{4}}{360}+2\cdot \frac{\pi^{4}}{72}=\frac{17\pi^{4}}{360}$

EDIT:

If I may expand somewhat on this sum using the same technique but a different f(z). Of course, it requires a couple known Euler sums as lemmata.

By considering $\displaystyle f(z)=\frac{(\gamma+\psi(-z))^{3}}{z^{2}}$, one can use the residues at 0 and the positive integers to find the sum.

Using the series for $\displaystyle(\gamma+\psi(-z))^{3}$ at z=n:

$\displaystyle \frac{1}{(z-n)^{3}}+\frac{3H_{n}}{(z-n)^{2}}+\frac{3(H_{n})^{2}}{z-n}-\frac{3H_{n}^{(2)}}{z-n}-\frac{\pi^{2}}{2(z-n)}+\cdot\cdot\cdot$

Thus, the residues at z=n are:

$\displaystyle\lim_{z\to n}\left(Res\left[\frac{1}{(z-n)^{3}}\cdot \frac{1}{z^{2}}\right]+Res\left[\frac{3H_{n}}{(z-n)^{2}}\cdot \frac{1}{z^{2}}\right]+Res\left[\frac{3(H_{n})^{2}}{z-n}\cdot \frac{1}{z^{2}}\right]-Res\left[\frac{H_{n}^{(2)}}{z-n}\cdot \frac{1}{z^{2}}\right]-Res\left[\frac{\pi^{2}}{2(z-n)}\right]\right)$

The first two require derivatives due to the pole at n being of order 3. But, we ultimately obtain the sums:

$\displaystyle 3\sum_{n=1}^{\infty}\frac{1}{n^{4}}-6\sum_{n=1}^{\infty}\frac{H_{n}}{n^{3}}+3\sum_{n=1}^{\infty}\frac{(H_{n})^{2}}{n^{2}}-3\sum_{n=1}^{\infty}\frac{H_{n}^{(2)}}{n^{2}}-\frac{\pi^{2}}{2}\sum_{n=1}^{\infty}\frac{1}{n^{2}}+\frac{\pi^{4}}{20}=0$

Also, the residue at z=0 is $\displaystyle\frac{\pi^{4}}{20}$, which can be found by using its Laurent expansion:

$\displaystyle f(z)=\frac{1}{z^{5}}-\frac{3\zeta(2)}{z^{3}}-\frac{3\zeta(3)}{z^{2}}+\frac{\pi^{4}}{20z}+\cdot\cdot\cdot$

Sum residues, evaluate known sums, call the quadratic sum being found H, set to 0 and solve for H.

$\displaystyle=\frac{\pi^{4}}{30}-\frac{\pi^{4}}{12}+3H-\frac{7\pi^{4}}{120}+\frac{\pi^{4}}{20}-\frac{\pi^{4}}{12}=0$

$\displaystyle \sum_{n=1}^{\infty}\frac{(H_{n})^{2}}{n^{2}}=\frac{17\pi^{4}}{360}$

Random Variable is an expert in this method and has refined it very well.

• I'm really sorry to bother you about such an old post, but may I ask how one can derive the identity $[1]$? Thank you very much. – M.N.C.E. Oct 28 '14 at 6:51
• One way is to use the generating functions and note that: $$1/3\frac{\log^{3}(1-x)}{x}dx+2\frac{Li_{3}(x)}{x}dx=2\sum_{n=1}^{\infty}\frac{H_{n}}{n^{2}}x^{n-1}+\sum_{n=1}^{\infty}\frac{H_{n}^{(2)}}{n}x^{n-1}-\sum_{n=1}^{\infty}\frac{H_{n}^{2}}{n}x^{n-1}$$. Then, integrate over 0 to t. – Cody Oct 28 '14 at 20:24
• I just used wolfram to integrate the left side to get $$1/3\log^{3}(1-t)\log(t)+\log^{2}(1-t)Li_{2}(1-t)-2\log(1-t)Li_{3}(1-t)+2Li_{4}(1-t)+2Li_{4}(t)-2\zeta(4)$$. Let $t\to 1$ and we get $$2\sum_{n=1}^{\infty}\frac{H_{n}}{n^{3}}+\sum_{n=1}^{\infty}\frac{H_{n}^{(2)}}{n^{2}}-\sum_{n=1}^{\infty}\frac{H_{n}^{2}}{n^{2}}=0$$. – Cody Oct 28 '14 at 20:34
• I see, thank you for clearing that up for me. – M.N.C.E. Oct 31 '14 at 2:18

Compute the generating function of the harmonic numbers: \begin{align} \sum_{n=1}^\infty H_nx^n &=\sum_{n=1}^\infty\sum_{k=1}^n\frac{x^n}{k}\\ &=\sum_{k=1}^\infty\sum_{n=k}^\infty\frac{x^n}{k}\\ &=\sum_{k=1}^\infty\sum_{n=0}^\infty\frac{x^{n+k}}{k}\\ &=-\frac{\log(1-x)}{1-x}\tag{1} \end{align} Integrating $(1)$ yields $$\sum_{n=1}^\infty\frac{H_{n-1}}{n}x^n=\frac12\log(1-x)^2\tag{2}$$ Therefore, $$\sum_{n=1}^\infty\frac{H_{n-1}}{n}e^{\pm2\pi inx}=\frac12\log(1-e^{\pm2\pi ix})^2\tag{3}$$ Multiplying and integrating gives \begin{align} \sum_{n=1}^\infty\frac{H_{n-1}^2}{n^2} &=\frac14\int_0^1\log(1-e^{2\pi ix})^2\log(1-e^{-2\pi ix})^2\,\mathrm{d}x\tag{4a}\\ &=\frac1{8\pi i}\oint\log(1-z)^2\log(1-1/z)^2\frac{\mathrm{d}z}{z}\tag{4b}\\ &=\frac1{8\pi i}\int_0^1\log(1-z)^2\left[-\pi i+\log(1-z)-\log(z)\right]^2\frac{\mathrm{d}z}{z}\\ &-\frac1{8\pi i}\int_0^1\log(1-z)^2\left[\pi i+\log(1-z)-\log(z)\right]^2\frac{\mathrm{d}z}{z}\tag{4c}\\ &=-\frac12\int_0^1\log(1-z)^2\left[\log(1-z)-\log(z)\right]\frac{\mathrm{d}z}{z}\tag{4d} \end{align} Explanation
$\mathrm{(4a)}$: multiply the conjugates of $(3)$ and integrate
$\mathrm{(4b)}$: convert to contour integral with $z=e^{2\pi ix}$
$\mathrm{(4c)}$: deflate the contour to lines above and below $[0,1]$
$\mathrm{(4d)}$: algebra

Contour $\color{#00A000}{\text{before}}$ and $\color{#C00000}{\text{after}}$ $\mathrm{(4c)}$:

$\hspace{4cm}$

Using $\log(1-z)=-u$, we get \begin{align} \int_0^1\log(1-z)^3\frac{\mathrm{d}z}{z} &=-\int_0^\infty u^3\frac{\mathrm{d}u}{e^u-1}\\ &=-\Gamma(4)\zeta(4)\\ &=-\frac{\pi^4}{15}\tag{5} \end{align} Using $\log(z)=-u$ and , we get \begin{align} \int_0^1\log(1-z)^2\log(z)\frac{\mathrm{d}z}{z} &=-\int_0^\infty\log(1-e^{-u})^2u\,\mathrm{d}u\tag{6a}\\ &=-2\sum_{n=1}^\infty\int_0^\infty\frac{H_{n-1}}{n}e^{-nu}u\,\mathrm{d}u\tag{6b}\\ &=-2\sum_{n=1}^\infty\frac{H_{n-1}}{n^3}\tag{6c}\\ &=\zeta(2)^2-3\zeta(4)\tag{6d}\\ &=-\frac{\pi^4}{180}\tag{6e} \end{align} Explanation
$\mathrm{(6a)}$: substitute $z=e^{-u}$
$\mathrm{(6b)}$: apply $(2)$
$\mathrm{(6c)}$: integrate
$\mathrm{(6d)}$: use this answer
$\mathrm{(6e)}$: evaluate

Combining $(4)$, $(5)$, and $(6)$ yields $$\sum_{n=1}^\infty\frac{H_{n-1}^2}{n^2}=\frac{11\pi^4}{360}\tag{7}$$ Noting that \begin{align} \sum_{n=1}^\infty\frac{H_{n-1}^2}{n^2} &=\sum_{n=1}^\infty\frac{\left(H_n-\frac1n\right)^2}{n^2}\\ &=\sum_{n=1}^\infty\left(\frac{H_n^2}{n^2}-2\frac{H_n}{n^3}+\frac1{n^4}\right)\tag{8} \end{align} we get, again using this answer, that \begin{align} \sum_{n=1}^\infty\frac{H_n^2}{n^2} &=\sum_{n=1}^\infty\frac{H_{n-1}^2}{n^2}+2\sum_{n=1}^\infty\frac{H_n}{n^3}-\zeta(4)\\ &=\frac{11\pi^4}{360}+5\zeta(4)-\zeta(2)^2-\zeta(4)\\ &=\frac{17\pi^4}{360}\tag{9} \end{align}

• This looks simpler than the other answers (may be a bit lengthy, but like solutions which don't require too much heavy machinery). +1 – Paramanand Singh Aug 25 '14 at 10:37


$$\mbox{Note that}\quad H_{n}=\int_{0}^{1}{1 - t^{n} \over 1 - t}\,\dd t =-n\int_{0}^{1}\ln\pars{1 - t}t^{n - 1}\,\dd t$$

Then, \begin{align} &\color{#c00000}{\sum_{n = 1}^{\infty}\pars{H_{n} \over n}^{2}} =\sum_{n = 1}^{\infty}\bracks{\int_{0}^{1}\ln\pars{1 - x}x^{n - 1}\,\dd x} \bracks{\int_{0}^{1}\ln\pars{1 - y}x^{n - 1}\,\dd y} \\[3mm]&=\int_{0}^{1}\int_{0}^{1} \ln\pars{1 - x}\ln\pars{1 - y}\sum_{n =1}^{\infty}\pars{xy}^{n - 1}\,\dd y\,\dd x \\[3mm]&=\int_{0}^{1}\ln\pars{1 - x} \color{#00f}{\int_{0}^{1}{\ln\pars{1 - y} \over 1 - xy}\,\dd y}\,\dd x\tag{1} \end{align}

\begin{align}&\color{#00f}{\int_{0}^{1}{\ln\pars{1 - y} \over 1 - xy}\,\dd y} =\int_{0}^{1}{\ln\pars{y} \over 1 - x\pars{1 - y}}\,\dd y =\int_{0}^{1}{\ln\pars{y} \over 1 - x + xy}\,\dd y \\[3mm]&=-\,{1 \over x}\int_{0}^{1}{\ln\pars{y} \over 1 - xy/\pars{x - 1}}\,{x\,\dd y \over x - 1} =-\,{1 \over x}\int_{0}^{x/\pars{x - 1}} {\ln\pars{\bracks{x - 1}y/x} \over 1 - y}\,\dd y \\[3mm]&=-\,{1 \over x}\int_{0}^{x/\pars{x - 1}}{\ln\pars{1 - y} \over y}\,\dd y ={1 \over x}\int_{0}^{x/\pars{x - 1}}{{\rm Li}_{1}\pars{y} \over y}\,\dd y \end{align} where $\ds{{\rm Li_{s}}\pars{z}}$ is the PolyLogarithm Function and $\ds{{\rm Li_{1}}\pars{z} = -\ln\pars{1 - z}}$. Hereafter, we'll use well known properties of them as reported in the above cited link: \begin{align}&\color{#00f}{\int_{0}^{1}{\ln\pars{1 - y} \over 1 - xy}\,\dd y} ={1 \over x}\int_{0}^{x/\pars{x - 1}}{\rm Li}_{2}'\pars{y}\,\dd y ={1 \over x}\,{\rm Li}_{2}\pars{x \over x - 1} \end{align}

Replacing the last result in expression $\pars{1}$: \begin{align} &\color{#c00000}{\sum_{n = 1}^{\infty}\pars{H_{n} \over n}^{2}} =\int_{0}^{1}\ln\pars{1 - x}\,{1 \over x}\,{\rm Li}_{2}\pars{x \over x - 1}\,\dd x =-\int_{0}^{1}{\rm Li}_{2}'\pars{x}{\rm Li}_{2}\pars{x \over x - 1}\,\dd x \\[3mm]&=-\int_{0}^{1}{\rm Li}_{2}'\pars{1 - x} {\rm Li}_{2}\pars{1 - {1 \over x}}\,\dd x =-\int_{0}^{1}{\rm Li}_{2}'\pars{1 - x}\bracks{-{\rm Li}_{2}\pars{1 - x} -\half\,\ln^{2}\pars{x}}\,\dd x \end{align} where we used Landen Identity. \begin{align} &\color{#c00000}{\sum_{n = 1}^{\infty}\pars{H_{n} \over n}^{2}} =\half\,{\rm Li}_{2}^{2}\pars{1} +\half\int_{0}^{1}{\rm Li}_{2}'\pars{1 - x}\ln^{2}\pars{x}\,\dd x \\[3mm]&={\pi^{4} \over 72} -\half\color{#00f}{\int_{0}^{1}{\ln^{3}\pars{x} \over 1 - x}\,\dd x} \quad\mbox{since}\quad{\rm Li}_{2}\pars{1} = {\pi^{2} \over 6}\tag{2} \end{align}

Finally, we have to evaluate the integral \begin{align}&\color{#00f}{\int_{0}^{1}{\ln^{3}\pars{x} \over 1 - x}\,\dd x} =\int_{0}^{1}\ln\pars{1 - x}\,\bracks{3\ln^{2}\pars{x}\,{1 \over x}}\,\dd x =-3\int_{0}^{1}{\rm Li}_{2}'\pars{x}\ln^{2}\pars{x}\,\dd x \\[3mm]&=3\int_{0}^{1}{\rm Li}_{2}\pars{x}\bracks{2\ln\pars{x}\,{1 \over x}}\,\dd x =6\int_{0}^{1}{\rm Li}_{3}'\pars{x}\ln\pars{x}\,\dd x \\[3mm]&=-6\int_{0}^{1}{\rm Li}_{3}\pars{x}\,{1 \over x}\,\dd x =-6\int_{0}^{1}{\rm Li}_{4}'\pars{x}\,\dd x=-6{\rm Li}_{4}\pars{1} =-6\zeta\pars{4}=-6\,{\pi^{4} \over 90}=\color{#00f}{-\,{\pi^{4} \over 15}} \end{align}

Replacing in $\pars{2}$: \begin{align} &\color{#66f}{\large\sum_{n = 1}^{\infty}\pars{H_{n} \over n}^{2}} ={\pi^{4} \over 72} - \half\,\pars{-\,{\pi^{4} \over 15}} =\color{#66f}{\large{17 \over 360}\,\pi^{4}} \end{align}

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.

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*}

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)$$

Different approach:

$$\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)$$

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}$$