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

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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

Answer 4

Using the divergent series $$ \sum_{j=0}^\infty(-1)^j\binom{2j}{n} =\mathrm{Re}\left[\frac{i^n}{(1-i)^{n+1}}\right]\tag{1} $$ which is gotten by computing $$ \lim_{\alpha\to1^-}\sum_{j=0}^\infty(-1)^j\binom{2j}{n}\alpha^{2j}\tag{2} $$ we get $$ \begin{align}\hspace{-1cm} &\sum_{k=1}^n\frac1{n^2+k^2}\\ &=\frac1{n^2}\sum_{k=1}^n\frac1{1+(k/n)^2}\\ &=\frac1{n^2}\sum_{k=1}^n\sum_{j=0}^\infty(-1)^j\frac{k^{2j}}{n^{2j}}\\ &=\sum_{j=0}^\infty\frac{(-1)^j}{n^{2j+2}}\sum_{k=1}^nk^{2j}\\ &=\sum_{j=0}^\infty\frac{(-1)^j}{n^{2j+2}}\left[\frac{n^{2j+1}}{2j+1}+\binom{2j}{0}\frac{n^{2j}}{2}+\binom{2j}{1}\frac{n^{2j-1}}{12}-\binom{2j}{3}\frac{n^{2j-3}}{120}+\dots\right]\\ &=\frac\pi{4n}-\frac1{4n^2}-\frac1{24n^3}+\frac1{2016n^7}-\frac1{4224n^{11}}+\frac1{1536n^{15}}+O\left(\frac1{n^{19}}\right)\tag{3} \end{align} $$ which is using $(1)$ to evaluate $$ \lim_{\alpha\to1^-}\sum_{k=1}^n\frac1{n^2+\alpha^2k^2}\tag{4} $$ Therefore $$ \begin{align} &\sum_{j,k=1}^n\frac1{j^2+k^2}-\sum_{j,k=1}^{n-1}\frac1{j^2+k^2}\\ &=\frac\pi{2n}-\frac1{n^2}-\frac1{12n^3}+\frac1{1008n^7}-\frac1{2112n^{11}}+\frac1{768n^{15}}+O\left(\frac1{n^{19}}\right)\tag{5} \end{align} $$ Applying the Euler-Maclaurin Sum Formula to $(5)$ we get $$ \begin{align} &\sum_{j,k=1}^n\frac1{j^2+k^2}\\ &=\frac\pi2\log(n)+C+\frac{4+\pi}{4n}-\frac{11+\pi}{24n^2}+\frac1{8n^3}+\frac{5+\pi}{240n^4}-\frac1{30n^5}-\frac{43+12\pi}{6048n^6}\\ &+\frac7{288n^7}+\frac{55+18\pi}{8640n^8}-\frac1{30n^9}-\frac{149+48\pi}{12672n^{10}}+\frac{29}{384n^{11}}+\frac{104845+33168\pi}{3144960n^{12}}\\ &-\frac{691}{2730n^{13}}-\frac{1411+448\pi}{10752n^{14}}+\frac{1793}{1536n^{15}}+\frac{1634635+520848\pi}{2350080n^{16}}-\frac{3617}{510n^{17}}\\ &+O\left(\frac1{n^{18}}\right)\tag{6} \end{align} $$ Computing the double sum for $n=10000$ and using $(6)$ we get that $C$ is approximately $$ -0.82586 11759 78831 08201 02008 35613 80953 63017 94512 34066 96955 08772 47540 16332 $$ which matches Kirill's answer.

Answer 5

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