Asymptotic estimate of double summation I would like to determine an asymptotic expansion for the following double summation:
$$\sum_{a=1}^{N/\sqrt {j}} \sum_{b=a}^{ja} \frac{1}{ab}$$
where $j$ is a real number $\geq 1$ and $N$ tends to $\infty$. In practice, the summation includes all pairs of integers  $a,b$ (with $a \leq b \leq ja $) such that the product $ab$ is $\leq N^2 $.
For $j=1$, the summation yields the infinite harmonic sum of squares, i.e. $\pi^2/6$. For $j>1$, the asymptotic expansion has the form $\log(N) \log(j) + O(1)$, where the $O(1)$ term converges to a value $k$ that depends on $j$. Interestingly, plotting $k$ vs $j$, the resulting function is discontinuous, with the most evident discontinuities occurring for $j$ integer. Is there any way to express this value explicitly? I am not necessarily searching a closed form, which may probably not exist. Rather, I would be interested to know potential alternative ways (e.g., using functions, series, and so on)  to express the value of $k$, different from the trivial definition given by the difference between the double summation and the log term. 
This question is a general formulation of a problem that is related to this and this other questions previously posted on MSE, for which no conclusive answer was provided despite multiple bounties. 
 A: Using the decomposition
$$
H(x) := \sum_{a \leq x} \frac{1}{a} = \log(x) + \gamma - \frac{\psi(x)}{x} + \int_{x}^{+\infty} \frac{\psi(t) d   t}{t^2},
$$
where $\psi(t) = \{ t \} - \frac{1}{2}$, one gets
$$
k(j) = - \sum_{n \geq 1} \frac{\psi(nj)}{n^2j } + \tilde{k}(j),
$$
where $\tilde{k}$ is an explicit Lipschitz function (see below). In particular :


*

*if $j$ is irrational, then $k$ is continuous at $j$.

*if $j = \frac{a}{q}$ is rational (with $a$ and $q$ coprime), then $k$ has a jump discontinuity, with $k(j^+) = k(j)$, and $k(j^+) - k(j^-) = \frac{\zeta(2)}{aq}$.


Moreover, $\tilde{k}$ has the form
$$
\tilde{k}(j)= -\frac{1}{2} \log(j)^2 + \gamma \log(j) + \frac{1}{j} \int_{1}^{+ \infty} \frac{\psi(jt) H(t) d t}{t^2}  + c,
$$
for some constant $c$.
A: Maple 2016 finds it for concrete values of $j$. For example,
asympt(sum(sum(1/(a*b), b = a .. 3*a, parametric), a = 1 .. N/sqrt(3)), N, 2)

$$\ln  \left( 3 \right) \ln  \left( N \right) +\ln  \left( 3 \right)
\ln  \left( 1/3\,\sqrt {3} \right) +1/18\,{\pi}^{2}+$$ $$\sum _{a=1}^{
\infty }\,{\frac {\Psi \left( a+2/3 \right) +\Psi \left( a+1/3
 \right) -2\,\Psi \left( a \right) }{3a}}+\ln  \left( 3 \right) \gamma +
{\frac {\,\ln  \left( 3 \right) \sqrt {3}-1/3\,\sqrt {3}}{2N}}+O
 \left( {N}^{-2} \right),$$
where  $\Psi(x)$ is the digamma function (http://www.maplesoft.com/support/help/Maple/view.aspx?path=Psi&term=Psi) and $\gamma=0.57721...$ is the Euler-Mascheroni constant.
PS. Also the Mathematica command 
Series[Sum[Sum[1/(a*b), {b, a, 3 a}], {a, 1, N/Sqrt[3]}, Assumptions -> N > 0], {N, Infinity, 1}]

does that in a somewhat more complicated form.
