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)$. In the chatroom we discussed about the asymptotic of $\displaystyle \sum_{i = 1}^n \sum_{j = 1}^n \frac1{i^2+j^2}$, and if we think of the inverse tangent integral, it's easy to see that $\displaystyle \sum_{i = 1}^n \sum_{j = 1}^n \frac1{i^2+j^2}\approx \operatorname{Ti}_2(n)\approx \frac{\pi}{2} \log(n)$ where I used
the well-known relation $\displaystyle \operatorname{Ti}_2(x)-\operatorname{Ti}_2(1/x)=\frac{\pi}{2}\operatorname{sgn}(x) \log|x|$. At this point, @robjohn posed the following  limit
$$\lim_{n\to\infty} \left(\displaystyle \sum_{i = 1}^n \sum_{j = 1}^n \frac1{i^2+j^2}-\frac{\pi}{2} \log(n)\right)$$
that looks like a pretty tough limit. Using $\coth(z)$ one can see the limit is approximately $-\dfrac{\pi^2}{12}$,
but how about finding a way to precisely compute the limit? 
 A: Here is a way a physicist gets a quick estimation:  
First we note that   
$$\sum_{i = 1}^n \sum_{j = 1}^n \frac1{i^2+j^2}\sim\int_0^{n}\int_0^{n}\frac{dxdy}{x^2+y^2}\;(n \to\infty)$$  
Now, the double integral  $I(n)$, because the integration is elementary I write only the final result:  
$$I(n)=\arctan(n)\ln(n)-\int_1^{n}\frac{\ln x}{1+x^2}dx-\int_{\frac{1}{n}}^{1}\frac{\arctan x}{x}dx$$  
Because the integrals in the last expression are bounded when $n \to\infty$, then  
$$I(n)\rightarrow \frac{\pi}{2}\ln n$$
A: 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.
A: 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.
A: Here is the beginning of an asymptotic expansion. First note that
$$
\begin{align}
\arctan\left(\frac kn\right)-\arctan\left(\frac {k-1}n\right)
&=\arctan\left(\frac{\frac1n}{1+\frac kn\frac {k-1}n}\right)\\
&=\arctan\left(\frac{n}{n^2+k^2-k}\right)\tag{1}
\end{align}
$$
Then
$$
\begin{align}\hspace{-1cm}
\sum_{k=1}^n\frac1{n^2+k^2}-\frac\pi{4n}
&=\sum_{k=1}^n\left[\frac1{n^2+k^2}-\frac1n\arctan\left(\frac{n}{n^2+k^2-k}\right)\right]\\
&=\sum_{k=1}^n\left[\frac1{n^2+k^2}-\frac1{n^2+k^2-k}+\frac{n^2}{3(n^2+k^2-k)^3}-\dots\right]\\
&=\sum_{k=1}^n\left[-\frac{k}{(n^2+k^2)(n^2+k^2-k)}+\frac{n^2}{3(n^2+k^2-k)^3}-\dots\right]\\
&\sim\int_0^1\left[-\frac1{n^2}\frac{x}{(1+x^2)^2}+\frac1{n^3}\frac1{3(1+x^2)^3}-\dots\right]\mathrm{d}x\\
&=-\frac1{4n^2}+O\left(\frac1{n^3}\right)\tag{2}
\end{align}
$$
Therefore,
$$
\sum_{k=1}^n\frac1{n^2+k^2}=\frac\pi{4n}-\frac1{4n^2}+O\left(\frac1{n^3}\right)\tag{3}
$$
By computing the edges of the square in index space, we get
$$
\begin{align}
\sum_{j,k=1}^n\frac1{j^2+k^2}-\sum_{j,k=1}^{n-1}\frac1{j^2+k^2}
&=2\sum_{k=1}^n\frac1{n^2+k^2}-\frac1{2n^2}\\
&=\frac\pi{2n}-\frac1{n^2}+O\left(\frac1{n^3}\right)\tag{4}
\end{align}
$$
The Euler-Maclaurin Sum Formula applied to $(4)$ gives
$$
\sum_{j,k=1}^n\frac1{j^2+k^2}=\frac\pi2\log(n)+C+\frac{\pi+4}{4n}+O\left(\frac1{n^2}\right)\tag{5}
$$
for some constant $C$.
Equation $(5)$ supports my comment that computationally, 
$$
\sum_{j,k=1}^n\frac1{j^2+k^2}-\frac\pi2\log(n)\sim C+\frac{\pi+4}{4n}\tag{6}
$$
where $C\approx-0.82586118$ and $\frac{\pi+4}{4}\approx1.785398$.
