Computing $\lim\limits_{n\to\infty} \Big(\sum\limits_{i = 1}^n \sum\limits_{j = 1}^n \frac1{i^2+j^2}-\frac{\pi}{2} \log(n)\Big)$.

I happened to struggle on the same problem 3 years ago. Here's another approach.

Start from $$\sum_{n=1}^\infty\frac1{x^2+n^2}=-\frac1{2x^2}+\frac{\pi}{2x}+\frac{\pi}{x(e^{2\pi x}-1)}$$ From this we have \begin{array} 1\sum_{m\le N}\sum_{n\le N}\frac{1}{m^2+n^2}&=\sum_{m\le N}\left(-\frac{1}{2m^2}+\frac{\pi}{2m}+\frac{\pi}{m(e^{2\pi m}-1)}-\sum_{n>N}\frac{1}{m^2+n^2}\right)\\ &=-\frac12\zeta(2)+\frac\pi2H_N+\pi\sum_{m\le N}\frac1{m(e^{2\pi m}-1)}-\sum_{m\le N}\sum_{n>N}\frac1{m^2+n^2}+O\left(\frac1N\right)\\ \end{array} Now $$\lim_{N\rightarrow\infty}\sum_{m\le N}\sum_{n>N}\frac1{m^2+n^2}=\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$$ and $$\sum_{m=1}^\infty\frac1{m(e^{2\pi m}-1)}=\sum_{m=1}^\infty\frac1m\sum_{k=1}^\infty e^{-2\pi mk}=-\sum_{k=1}^\infty \log(1-e^{-2\pi k})=-\log(\eta(i)e^{\pi/12})$$ The value of $\eta(i)$ is known to be $\Gamma(1/4)/(2\pi^{3/4})$.

Hence we have \begin{array} 1&\lim_{N\rightarrow\infty}\left(\sum_{m\le N}\sum_{n\le N}\frac{1}{m^2+n^2}-\frac\pi2\log N\right)\\ &=\pi\left(-\log\Gamma\left(\frac14\right)-\frac\pi{12}+\log2+\frac34\log\pi\right)-\frac{\zeta(2)}2+\frac{\pi\gamma}2-G\\ &=\frac14\pi\log\left(\frac{16\pi^3e^{2\gamma}}{\Gamma(\frac14)^4}\right)-G-\zeta(2) \end{array} same as Kirill's answer.

The limit in question is equal to $$\def\tfrac#1#2{{\textstyle\frac{#1}{#2}}} \tfrac14\pi\log\left(\frac{16\pi^3e^{2\gamma}}{\Gamma(\frac14)^4}\right) -G-\zeta(2) \\ = -0.82586\ 11759\ 78831\ 08201\ 02008\ 35613\ 80953\ 63017\ 94512\ 34066\ 96955\ 08772 \ldots $$ in terms of Euler's $\gamma$ and Catalan $G$.

Demonstration. The sum $$ S(N) = \sum_{1\leq m,n\leq N} \frac1{m^2+n^2} \sim \tfrac12\pi\log N $$ is in an inconvenient form, because it is taken over the square $[1,N]^2$ instead of the circle of radius $N$: $$ S_2(N) = \sum_{m,n\geq 1}\frac{[m^2+n^2\leq N^2]}{m^2+n^2} \sim \tfrac12\pi\log N. $$ This is so, because the circle-sum can be rewritten as a one-dimensional sum involving the function $r_2(k)$ that counts the number of ways to write $k$ as a sum of two squares of (signed) integers. So we can write the sum $S_2$ as a quarter of the sum over a circle in $\mathbb{Z}^2$ excluding the axes, like so: $$\begin{eqnarray} S_2(N) &=& \frac14 \sum_{m,n\in\mathbb{Z}}' \frac{[m^2+n^2\leq N^2]}{m^2+n^2} - H_N^{(2)} \\&=& \frac14\sum_{1\leq n\leq N^2} \frac{r_2(n)}{n} - H_N^{(2)}. \end{eqnarray} $$ Here $H_N^{(2)}$ has the limit $\zeta(2)=\frac{\pi^2}{6}$, and the prime above the sum means the sum omits $(m,n)=(0,0)$.

The sum $$ g(n) = \sum_{1\leq k\leq n} \frac{r_2(k)}{k} = \pi\log N + \pi S + O(N^{-1/2}) $$ has a standard asymptotic form in terms of the Sierpinski constant $$ S = \log \left(\frac{4\pi^3e^{2\gamma}}{\Gamma(\frac14)^4}\right). $$

The difference $S(N)-S_2(N)$ between the circle-sum and the square-sum does not diverge as $N\to\infty$, and we may sandwich it above and below with two integrals that converge to the same value of $$ \int_1^N dx\int_1^N dy \frac{[x^2+y^2>N^2]}{x^2+y^2}. $$ This integral can be computed in polar coordinates, by writing $$ \begin{eqnarray} 2\int_{N}^{N\sqrt{2}}\frac{dr}{r} \int_0^{\pi/4}d\phi\,[r\cos\phi < N] &=& 2\int_1^{\sqrt{2}}\frac{dr}{r}\left(\frac\pi4 - \arccos\frac1r\right) \\ &=& \tfrac12\pi\log2 - G. \end{eqnarray}$$

Putting everything together leads to the expression given above.