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)$.

• Presumably the size of the discontinuity at any rational number $J=c/d$ in lowest terms is equal to the contribution from the infinitely many points on the line $y=Jx$, which is $\pi^2/6cd$. Presumably also this function of $J$ is continuous at irrational values of $J$. But that's not a particularly helpful statement.... – Greg Martin May 29 '14 at 22:25
• Can you share your plot with us? – cactus314 Jun 2 '14 at 10:25
• Anatoly - How sure are you about your asymptotic expansion? For $J, K$ positive integers, I'm getting $\log(J)\log(K) + \log^2(K)/2 + O(1+\log(JK))$. – John M Jun 2 '14 at 19:53
• This thing reminds me on the 1D-quasicrystal case. Some 'tricks' working on this would include constructing the lattice using a product of (a) the 2D version of the Dirac comb function (b) a z-rotated Heavyside stripe and (c) $1/xy$. That can be easily Fourier-transformed. – Rudi_Birnbaum Aug 9 '16 at 11:44
• This may be an intriguing similarity. Could you elaborate on it? – Anatoly Aug 9 '16 at 12:53

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).
• Thank you for your nice answer. The specific cases where $J$ is integer are actually solvable in several ways. For example, we could use the formulas previously reported in an interesting paper by O. Kouba (see arxiv.org/pdf/1010.1842.pdf) to get closed forms. However, the difficulty in this problem is to solve it for non-integer values of $J$ (and particularly for irrational values). – Anatoly Oct 9 '14 at 13:57
• That paper gives a formula for another sum, namely $\displaystyle \sum_{n=1}^\infty \frac{\log(J)-(H_{Jn}-H_n)}{n}$ with $J$ integer. Making some calculations, we can use this result to get closed forms for the constant term of our double harmonic summation for $J$ integer. In particular, the case $J=2$ has the expansion $\displaystyle \log(K)\log(2)+ \frac{\pi^2}{12} + \gamma \log(2) + \log^2(2)$. The sum of the 3 last terms is the constant term, equal to $1.703015...$, in accordance with the Figure. Unfortunately, it seems that the method of Kouba cannot be applied to non-integer $J$. – Anatoly Oct 10 '14 at 0:05