Double harmonic summation

Let us consider a lattice formed by all points with integer and positive coordinates on a Cartesian plane, and where $K$ is the maximal value for the $x$-axis. Let us assign to each lattice point the value $1/(xy)$, where $x$ and $y$ are its coordinates. Drawing the two lines $y = x$ and $y = Jx$ (with $J$ real number $>1$), a triangular region between the two lines is identified. I would like to determine an asymptotic expansion for the summation of the values of all points included in this triangular region when $K$ tends to infinity.

I calculated that this asymptotic expansion is given by $\log(K) \log(J) + O(1)$. However, I am particularly interested in the limit of this constant term when $K$ tends to infinity, since it shows an irregular behaviour with several discontinuities when plotted against $J$ (see the Figure below, showing the value of this constant term on the $y$-axis vs $J$ on the $x$-axis).

As shown in the Figure, for $J = 1$ (where the two lines coincide and the triangle reduces to a line crossing all points with $x = y$), the limit of the constant term, and in this case of the whole summmation, corresponds to the infinite sum of inverse squares, yielding $\pi^{2}/6 \approx 1.64$. Starting from this value, the function then shows discontinuities over its whole range. Just to describe two of them that are evident in magnitude and that have been highlighted in the Figure: for $J = 1.49999$ the limit of the constant term is about $1.025$, but for $J = 1.5$ it abruptly increases to about $1.30$. Similarly, for $J = 1.99999$ the limit of the constant term decreases to about $0.88$, whereas for $J = 2$ it abruptly increases to about $1.70$. The reason for all these discontinuities is that, as the $y = Jx$ line rotates counterclockwise (i.e., as $J$ increases), new sets of aligned points are progressively and discontinuously included in the triangular region. Accordingly, the largest discontinuities occur for integer values of $J$, as for these values new and relatively large sets of points are included in the region. Generalizing, the magnitude of each jump is equal to the sum of all points corresponding to the slope $J=\frac{y}{x}$ (where $\frac{y}{x}$ is an irreducible fraction), which in turn is given by the sum of inverse squares multiplied by $\displaystyle \frac{1}{xy}$. In this view, this question is equivalent to that of determining an asymptotic expansion for the double harmonic sum $\displaystyle \frac{1}{xy}$, calculated over all possibile coprime values of $x$ and $y$ such that $\frac{y}{x}$ is a rational number included in the interval between $1$ and $J$.

I would like to determine a general expression for the constant term of this summation. I tried a number of approaches, including application of Euler-MacLaurin formula or Fourier transformations, but each attempt failed. I am not necessarily searching a closed form in the strict sense, which probably could not exist. Rather, I would be very interested in finding a more compact or elegant general expression for this term, in comparison with the trivial way of reporting it as the difference between the double harmonic summation and the logarithmic term $\log(K) \log(J)$.

Solution 1:

Just a partial answer for now. We have: $$\sum_{\substack{x\leq y\leq Jx\\1\leq x\leq N}}\frac{1}{xy}=\sum_{x=1}^{N}\sum_{y=x}^{\lfloor Jx\rfloor}\frac{1}{xy}=\sum_{x=1}^{N}\frac{H_{\lfloor Jx\rfloor}-H_{x-1}}{x}$$ and, through summation by parts: $$\sum_{n=1}^{N}\frac{H_{n-1}}{n}=-H_{N}^{(2)}+\sum_{n=1}^{N}\frac{H_n}{n}=-H_N^{(2)}+H_N^2-\sum_{n=1}^{N-1}\frac{H_n}{n+1}$$ so: $$\sum_{n=1}^{N}\frac{H_{n-1}}{n}=\sum_{1\leq i<j\leq N}\frac{1}{ij}=\frac{1}{2}\left(H_N^2-H_N^{(2)}\right).\tag{1}$$ On the other hand, the identity $H_n=\gamma+\Psi(n+1)$, where $\Psi$ is the digamma function, gives: $$\sum_{x=1}^{N}\frac{H_{\lfloor Jx\rfloor}}{x}=\gamma H_n+\sum_{x=1}^{N}\frac{\Psi(\lfloor Jx\rfloor+1)}{x}.\tag{2}$$ In the special case $J=2$ we have: $$\begin{eqnarray*}\sum_{n=1}^{N}\frac{H_{2n}}{n}&=&H_N H_{2N}-\sum_{n=1}^{N-1}H_n\left(\frac{1}{2n+2}+\frac{1}{2n+1}\right)\\&=&H_N H_{2N}-\frac{1}{4}\left(H_N^2-H_N^{(2)}\right)-\sum_{n=1}^{N-1}\frac{H_n}{2n+1}\\&=&H_N H_{2N}-\frac{1}{2}\left(H_N^2-H_N^{(2)}\right)-\sum_{n=1}^{N-1}\frac{H_n}{(2n+1)(2n+2)}\\&=&H_N\left(H_{2N}-\frac{1}{2}H_N\right)+\log^2 2+O\left(\frac{\log N}{N^2}\right)\\&=&\frac{1}{2}\log^2 N + (\log 2+\gamma)\log N+\left(\log^2 2+\gamma\log 2+\frac{\gamma^2}{2}\right)+O\left(\frac{\log N}{N}\right)\tag{3}\end{eqnarray*}$$ and for integer values of $J$ the asymptotics of $\sum_{n=1}^{N}\frac{H_{Jn}}{n}$ can be computed with the same technique (summation by parts).