<s>[INST] On the asymptotic behaviour of the Cauchy product of harmonic numbers $\sum_{k=1}^n H_k H_{n-k+1}$ In this post we take in our hands two simple tools. The first is the generating function of the harmonic numbers, you can see it in this Wikipedia, section 3, that holds for $|z|<1$. The second is the Cauchy product formula $$ \left( \sum_{n=1}^\infty a_n \right)  \left( \sum_{n=1}^\infty b_n \right) =  \sum_{n=1}^\infty a_k b_{n-k+1},$$ where the convergence is assumed, there are theorems that tell you when is convergent.
In our case our factors are convergents because are the same, and this is well defined: we consider the square of the generating function and after we take the integral $\int_0^{1/2}$ to get if there are no mistakes $$\sum_{n=1}^\infty\frac{1}{(n+2)2^{n+2}}\sum_{k=1}^{n}H_k H_{n-k+1}=\int_0^{1/2}\left(\sum_{n=1}^\infty H_n z^n\right)^2 dz=\int_0^{1/2}\left(\frac{\log(1-z)}{1-z}\right)^2 dz.$$ Notice that is required use the Cauchy product and swap the sign of the series and integral. The integral in RHS is computed in a closed form as $2(\log 2-1)^2$. 
I don't know if is well known from the literature this sequence $$\sum_{k=1}^n H_k H_{n-k+1}.$$ The sequence starts as $1,3, \frac{71}{12}, \frac{29}{3}, \frac{638}{45}, \frac{349}{18}, \frac{14139}{560}, \frac{79913}{2520}\ldots$ And I would like to know what's about its asymptotic behaviour. One know that the plot of this arithmetic function is smooth (with a well defined slope), see how  Wolfram Alpha can show us the plot of partial sums if you type the code 
sum HarmonicNumber[k]HarmonicNumber[1000-k+1], from k=1 to 1000
in its online calculator.

Question. What's about the asymptotic behaviour of $$\sum_{k=1}^{n}H_k H_{n-k+1}$$
  as $n\to\infty$? I am saying a big oh, or small oh statement or your answer as an asymptotic equivalence. You can provide me hints to get it with summation (I don't know if it is easy), or have you another idea? Thanks in advance. 

Thus feel free to add hints, references if you need it, or a more detailed answer.
 [/INST] 
Lemma For $n \ge 2$,
  $$\begin{align}
\rho_n \stackrel{def}{=}\sum_{k=1}^{n-1} H_k H_{n-k} 
&= (n+1)\left[\psi'(n+2)-\psi'(2) + (\psi(n+2)-\psi(2))^2 \right]\\
&= (n+1)\left[\psi'(n+2)-\left(\frac{\pi^2}{6}-1\right) + (H_{n+1} - 1)^2\right]
\end{align}
$$
  where $\psi(x)$ is the digamma function.

In following expansion
$$
\rho_n = \sum_{k=1}^{n-1} H_k H_{n-k}
= \sum_{k=1}^{n-1}\left(\sum\limits_{p=1}^k \frac{1}{p}\right)\left(\sum_{q=1}^{n-k}\frac{1}{q}\right)
$$
the term $\displaystyle\;\frac{1}{pq}$ appear $(n+1) - (p+q)$ times. This leads to
$$\begin{align}
\rho_n
&= (n+1)\sum_{\substack{1 \le p, q\\p+q \le n}}\frac{1}{pq} - 2\sum_{\substack{1 \le p, q\\p+q \le n}}\frac{1}{p}
= (n+1)\sum_{\substack{1 \le p, q\\p+q \le n}}\frac{1}{pq} - 2\sum_{p=1}^{n-1}\frac{n-p}{p}\\
&= (n+1)\sum_{\substack{1 \le p, q\\p+q \le n}}\frac{1}{pq} - 2n(H_{n} - 1)\\
\implies\quad \frac{\rho_n}{n+1} - \frac{\rho_{n-1}}{n} &= 
\sum_{\substack{1 \le p, q\\p+q = n}}\frac{1}{pq} - \frac{2n}{n+1}(H_{n} - 1) + \frac{2(n-1)}{n}(H_{n-1}-1)\\
&= \frac{1}{n}\sum_{\substack{1 \le p, q\\p+q = n}}\left(\frac{1}{p}+\frac{1}{q}\right) + \frac{2}{n+1}(H_n-1)  + \frac{2}{n} H_{n-1} = \frac{2}{n+1}(H_{n} - 1)
\end{align}
$$
Notice last expression can be rewritten as
$$\frac{2}{n+1}(H_{n} - 1) = 2(H_{n+1}-H_n)(H_n-1) = (H_{n+1}-1)^2 - (H_n-1)^2  - \frac{1}{(n+1)^2}\\
= (\psi'(n+2) + (\psi(n+2)-\psi(2))^2)
- (\psi'(n+1) + (\psi(n+1)-\psi(2))^2)$$
The expression
$\displaystyle\;\frac{\rho_n}{n+1} - (\psi'(n+2) + (\psi(n+2)-\psi(2))^2)\;$
is a constant independent of $n$.
Evaluate this expression $n = 2$, we find the constant equal to $-\psi'(2)$ and our lemma follows.
For large $x$, we have following asymptotic expansion of $\psi$ and $\psi'$.
$$
\psi(x) \asymp 
\log(x) - \frac{1}{2x} - \sum_{k=1}^\infty \frac{B_{2k}}{2k\;x^{2k}}
\quad\text{ and }\quad
\psi'(x) \asymp \frac{1}{x} + \frac{1}{2x^2} + \sum_{k=1}^\infty \frac{B_{2k}}{x^{2k+1}}
$$
where $B_k$ is the $k^{th}$ Bernoulli number. As a result,
$$\rho_{n-2} \asymp (n-1)\left[\small
\frac{1}{n} + \frac{1}{2n^2} + \sum_{k=1}^\infty \frac{B_{2k}}{n^{2k+1}} -
\left(\frac{\pi^2}{6} - 1\right)
+  \left(
\log(n) - \frac{1}{2n} - \sum_{k=1}^\infty \frac{B_{2k}}{2k\;n^{2k}} - (1-\gamma)
\right)^2
\right]
$$
For large $n$, the leading behavior of $\rho_n$ is $n\log(n)^2$.
</s>