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.
 A: 
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$.
A: 
We show the  following  asymptotic expansion is valid
\begin{align*}
\sum_{k=1}^{n-1}H_kH_{n-k}\sim n\log^2 n+2(\gamma -1)n\log n+\left(\gamma ^2-2\gamma+2-\frac{\pi^2}{6}\right)n\tag{1}
\end{align*}
with $\gamma=0.57721\ldots$ the Euler-Mascheroni constant.

The answer is based upon Singularity Analysis which describes the asymptotic behavior of a large class of functions with moderate growth and decay and the strongly connected asymptotic behavior of the coefficients of their generating functions.
The functions under consideration have as central elements of an expansion at a singularity $\zeta$ 
\begin{align*}
\left(1-\frac{z}{\zeta}\right)^{-\alpha}\left(\log\frac{1}{1-\frac{z}{\zeta}}\right)^\beta\qquad\qquad \alpha,\beta\in\mathbb{C}
\end{align*}
These elements contribute asymptotically  terms of type
\begin{align*}
\zeta^{-n}n^{\alpha-1}\left(\log n\right)^\beta
\end{align*}
and the generating function of the sum in (1) is precisely of this type.
Note: The basic facts mentioned here are from chap. VI of  the classic Analytic Combinatorics  by P. Flajolet and R. Sedgewick.

Since 
  \begin{align*}
\sum_{n=1}^\infty H_n z^n=\frac{1}{1-z}\log\frac{1}{1-z}
\end{align*}
  we obtain as generating function of the sum   in  (1)
  \begin{align*}
(1-z)^{-2}\left(\log \frac{1}{1-z}\right)^2&=\left(\sum_{n=1}^\infty     H_n z^n\right)^2
=\sum_{n=2}^\infty\left(\sum_{{k+l=n}\atop{k,l\geq 1}} H_kH_l\right)z^n\\
&=\sum_{n=2}^\infty\left(\sum_{k=1}^{n-1} H_kH_{n-k}\right)z^n\tag{2}
\end{align*}

We can now apply the following theorem from the book:

Theorem VI.2: Let $\alpha$ be an arbitrary complex number in $\mathbb{C}\setminus\mathbb{Z}_{\leq 0}$. The coefficient of $z^n$ in the    function
  \begin{align*}
f(z)=(1-z)^{-\alpha}\left(\frac{1}{z}\log\frac{1}{1-z}\right)^{\beta}\tag{3}
\end{align*}
  admits for large $n$ a full asymptotic expansion in descending powers of $\log n$,
  \begin{align*}
f_n\equiv[z^n]f(z)\sim\frac{n^{\alpha-1}}{\Gamma{(\alpha)}}\left(\log n\right)^\beta
\left[1+\frac{C_1}{\log n}+\frac{C_2}{\log^2  n}+\cdots\right]\tag{4},
\end{align*}
  where 
  \begin{align*}
\left.C_k=\binom{\beta}{k}\Gamma(\alpha)\frac{d^k}{ds^k}\frac{1}{\Gamma(s)}\right|_{s=\alpha}
\end{align*}
Note: A coefficient of $\frac{1}{z}$ is introduced in front of the logarithm in (3) since $\log(1-z)^{-1}=z+O(z^2)$: in this  way, $f(z)$ is a bona fide powers series in $z$, even when $\beta$ is not an integer. Such a factor does not affect the asymptotic expansions in a logarithmic scale near $z=1$.

$$ $$

In order to derive the asymptotic expansion of the sum in (1) we apply (4) with $\alpha=2,\beta=2$ and obtain
  \begin{align*}
\sum_{k=1}^{n-1} H_kH_{n-k}&=[z^n](1-z)^{-2}\left(\log \frac{1}{1-z}\right)^2\\
&\sim n\left(\log n\right)^2\left(1+\frac{C_1}{\log  n}+\frac{C_2}{\log^2  n}+\cdots\right)\tag{5}
\end{align*}
Using the Digamma function
  \begin{align*}
\psi(s)=\frac{\Gamma^\prime(s)}{\Gamma(s)}\qquad s\in \mathbb{C}\qquad\text{and}\qquad
\psi(n)=H_n-\gamma\qquad n\in\mathbb{N}
\end{align*}
  we obtain from
  \begin{align*}
\left.C_k=\binom{2}{k}\frac{d^k}{ds^k}\frac{1}{\Gamma(s)}\right|_{s=2}
\end{align*}
  and some help of Wolfram Alpha 
  \begin{align*}
C_1&=\left.\binom{2}{1}\frac{d}{ds}\frac{1}{\Gamma(s)}\right|_{s=2}=\left.-2\frac{\psi(s)}{\Gamma(s)}\right|_{s=2}\
=\gamma-1\\
\\
C_2&=\left.\frac{d^2}{ds^2}\frac{1}{\Gamma(s)}\right|_{s=2}=\left.\frac{\left(\psi(s)\right)^2-\psi^\prime(s)}{\Gamma(s)}\right|_{s=2}
=\gamma^2-2\gamma+2-\frac{\pi^2}{6}
\end{align*}
  Putting $C_1$ and $C_2$ in (5) and the claim (1) follows.

This result can also be found in figure VI.5 on p.388. Note that  $C_k=0$ for $k>2$. So  we don't  get any more information from this theorem. 
