An expansion of double sum Let   $$s_n=\sum_{j=1}^n\sum_{i=n+1}^{+\infty} \dfrac{1}{i^2+j^2}$$
I need to find an expansion as $s_n=a+\frac bn+ o(\frac 1n)$. Numerically $a=G$ (Catalan's constant) and $b=-\frac 12$.
I can show that $a=G$:
By Comparing discrete sums and integrals  $$\lim_{n\rightarrow\infty}s_n =\lim_{n\rightarrow\infty}\int_n^{+\infty}\int_0^n\frac{\mathrm{d}x\,\mathrm{d}y}{x^2+y^2}=\int_1^\infty\frac1t\tan^{-1}(\frac1t)\,\mathrm{d}t=G$$
I need help to show that $$\lim_{n\rightarrow\infty} n(s_n-G)=-\frac 12$$
 A: Let
$$f(x,y) = \frac{1}{x^2 + y^2}.$$
The series $s_n$ is the Riemann sum approximation to the integral
$$G = \int_0^1 \int_{1}^\infty f(x,y) dydx.$$
with the sum
$$s_n = \sum_{i = 1}^n\sum_{j = n + 1}^{\infty} \frac{1}{n^2}f\left(\frac{i}{n}, \frac{j}{n}\right)$$
In each small square $[(i - 1) / n, i/n] \times [(j - 1 / n, j/n)]$, we have the first order approximation
$$\int_{(i - 1) / n}^{i / n} \int_{(j - 1) / n}^{j / n} f(x,y)dydx - \frac{1}{n^2} f\left(\frac{i}{n}, \frac{j}{n}\right)$$
$$= \int_{(i - 1) / n}^{i / n} \int_{(j - 1) / n}^{j / n} \partial_x f\left(\frac{i}{n}, \frac{j}{n}\right) \left(x - \frac{i}{n}\right) + \partial_y f\left(\frac{i}{n}, \frac{j}{n}\right) \left(y - \frac{i}{n}\right) dydx + O(n^{-4})$$
$$= -\frac{1}{2n^3}(\partial_x f\left(\frac{i}{n}, \frac{j}{n}\right) + \partial_y f\left(\frac{i}{n}, \frac{j}{n}\right)) + O(n^{-4}).$$
Thus we conclude that
$$\lim_{n \to \infty} n(G - s_n) = -\frac{1}{2} \int_0^1 \int_{1}^\infty \partial_x f(x, y) + \partial_y f(x, y) dydx.$$
By FTC, we can easily evaluate the integral to be $-1$, as
$$\int_0^1 \int_{1}^\infty \partial_x f(x, y) dydx = \int_{1}^\infty f(1, y) - f(0, y)dy = \frac{\pi}{4} - 1,$$
$$\int_0^1 \int_{1}^\infty \partial_y f(x, y) dydx = \int_{0}^1 -f(x, 1)dx = -\frac{\pi}{4}.$$
So we have the desired result.
