Infinite Series $\sum\limits_{n=1}^\infty\left(\frac{H_n}n\right)^2$

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.

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.

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