I'm looking to approximate the behavior of the following quantity, as a function of $n$ and $s$
$$f(s,n)=\frac{1}{n^2 (n + 2)}\sum_{i,j=1}^n (h_i-h_j)^2$$ where $$h_i=\left(1 - \frac{1}{i}\right)^{2 s}\frac{1}{i}$$
For $s=0$, this is well approximated by $\frac{\pi^2}{3n^2}$, derived here. How does this behave for $s>0$?
Background: This is the variance of loss after applying $s$ steps of gradient descent and learning rate 1 to minimize quadratic form with eigenvalues $1,\frac{1}{2},\frac{1}{3},\ldots$ starting with random initial point on the unit circle.


One approach is to write \begin{align}g(n,s)&=\frac12n^2(n+2)f(n,s)=\frac12\sum_{i,j=1}^n (h_i-h_j)^2\\&=n\sum_{k=1}^n\left(1-\frac1k\right)^{4s}\frac1{k^2}-\left(\sum_{k=1}^n\left(1-\frac1k\right)^{2s}\frac1k\right)^2\\&=n\sum_{r=0}^{4s}\binom{4s}r(-1)^rH_n^{(2+r)}-\left(\sum_{t=0}^{2s}\binom{2s}t(-1)^tH_n^{(1+t)}\right)^2\end{align} where $H_n^{(m)}$ is the generalised harmonic number. Use Euler-Maclaurin to derive the result $$H_n^{(m)}\sim\zeta(m)+\frac{n^{1-m}}{1-m}+\frac{n^{-m}}2-\sum_{p=1}^\infty\binom{m+2p-1}{2p}\frac{B_{2p}n^{-m-2p+1}}{m+2p-1}$$ from which we can easily cut off higher-order terms. Up to ${\cal O}(n^{-3})$, we have \begin{align}g(n,s)&\sim\small n\left(\sum_{r=0}^{4s}\binom{4s}r(-1)^r\zeta(2+r)-n^{-1}+\frac{n^{-2}}2-\frac{n^{-3}}6-4s\left(-\frac{n^{-2}}2+\frac{n^{-3}}2\right)+\frac{4s(4s-1)}2\left(-\frac{n^{-3}}3\right)\right)\\&\quad\small-\left(\log n+\gamma+\frac{n^{-1}}2-\frac{n^{-2}}{12}+\sum_{t=1}^{2s}\binom{2s}t(-1)^t\zeta(1+t)-2s\left(-n^{-1}+\frac{n^{-2}}2\right)+\frac{2s(2s-1)}2\left(-\frac{n^{-2}}2\right)\right)^2\\&\sim \small n\sum_{r=0}^{4s}\binom{4s}r(-1)^r\zeta(2+r)-1+\frac{(4s+1)n^{-1}}2-\frac{(4s+1)^2n^{-2}}6\\&\quad\small-\left(\log n+\gamma+\sum_{t=1}^{2s}\binom{2s}t(-1)^t\zeta(1+t)+\frac{(4s+1)n^{-1}}2-\frac{(12s^2-6s+1)n^{-2}}{12}\right)^2.\end{align} Thus \begin{align}g(n,s)&\sim a_1n-(1+a_2^2)+(b_1-2a_2b_2)n^{-1}+(c_1-2a_2c_2-b_2^2)n^{-2}\\&\quad-\log^2n-2(a_2+b_2n^{-1}+c_2n^{-2})\log n\end{align} where \begin{align}a_1=\sum_{r=0}^{4s}\binom{4s}r(-1)^r\zeta(2+r)&\quad b_1=\frac{4s+1}2&c_1=-\frac{(4s+1)^2}6\\a_2=\gamma+\sum_{t=1}^{2s}\binom{2s}t(-1)^t\zeta(1+t)&\quad b_2=\frac{4s+1}2&c_2=-\frac{12s^2-6s+1}{12}\end{align} so that $f(n,s)\sim2n^{-3}g(n,s)$.