$\sum_{i,j=1}^N \frac{1}{(i^2 + j^2)}$ (i..e. $1/r^2$ over the square lattice $i,j\in\{1,2,\ldots,N\}$) Does anyone have any insight as to how I might find a closed-form expression for
$\sum_{i=1}^N \sum_{j=1}^N\frac{1}{(i^2 + j^2)}$ ?
It feels like it should be straightforward because it's easy to do as an integral in polar coordinates if you include the axes. But here the summation is over a square lattice, and excludes the axis. Wolfram Alpha can't do it, and neither can I!

An asymptotic formula for large $N$ may be derived as follows. First,
$$
\sum\limits_{m = 1}^N {\frac{1}{{m^2  + n^2 }}}  =  - \frac{1}{{2n^2 }} + \frac{\pi }{{2n}}\coth (\pi n) - \frac{{\operatorname{Im} \psi (1 + N + \mathrm{i}n)}}{n}
$$
where $\psi$ is the digamma function. From the known asymptotics of the digamma function,
\begin{align*}
\operatorname{Im} \psi (1 + N + \mathrm{i}n) &= \arg (N + \mathrm{i}n) + \frac{n}{{2(n^2  + N^2 )}} + \mathcal{O}\!\left( {\frac{1}{{N^2 }}} \right) \\ & = \arctan \left( {\frac{n}{N}} \right) + \frac{n}{{2(n^2  + N^2 )}} + \mathcal{O}\!\left( {\frac{1}{{N^2 }}} \right).
\end{align*}
Now
\begin{align*}
\sum\limits_{n = 1}^N {\frac{1}{n}\arctan \left( {\frac{n}{N}} \right)} & = \frac{1}{N}\sum\limits_{n = 1}^N {\frac{N}{n}\arctan \left( {\frac{n}{N}} \right)} \\ & = \int_0^1 {\frac{{\arctan x}}{x}\mathrm{d}x}  + \mathcal{O}\!\left( {\frac{1}{N}} \right) = G + \mathcal{O}\!\left( {\frac{1}{N}} \right)
\end{align*}
where $G$ is Catalan's constant. Also
$$
\sum\limits_{n = 1}^N {\frac{1}{{2(n^2  + N^2 )}}}  = \mathcal{O}\!\left( {\frac{1}{N}} \right),\;\, \sum\limits_{n = 1}^N {\mathcal{O}\!\left( {\frac{1}{{N^2 }}} \right)}  = \mathcal{O}\!\left( {\frac{1}{N}} \right),\;\, \sum\limits_{n = 1}^N { - \frac{1}{{2n^2 }}}  =  - \frac{{\pi ^2 }}{{12}} + \mathcal{O}\!\left( {\frac{1}{N}} \right)
$$
and
\begin{align*}
\sum\limits_{n = 1}^N {\frac{\pi }{{2n}}\coth (\pi n)} & = \frac{\pi }{2}\sum\limits_{n = 1}^N {\frac{1}{n}}  + \frac{\pi }{2}\sum\limits_{n = 1}^N {\frac{{\coth (\pi n) - 1}}{n}} \\ & = \frac{\pi }{2}\log N + \frac{\pi }{2}\gamma  + \frac{\pi }{2}\sum\limits_{n = 1}^\infty  {\frac{{\coth (\pi n) - 1}}{n}}  + \mathcal{O}\!\left( {\frac{1}{N}} \right),
\end{align*}
using the asymptotics of the harmonic numbers ($\gamma$ is the Euler–Mascheroni constant). In summary
$$
\sum\limits_{n = 1}^N {\sum\limits_{m = 1}^N {\frac{1}{{m^2  + n^2 }}} }  = \frac{\pi }{2}\log N + C + \mathcal{O}\!\left( {\frac{1}{N}} \right)
$$
with
$$
C = \frac{\pi }{2}\gamma  - G - \frac{{\pi ^2 }}{{12}} + \frac{\pi }{2}\sum\limits_{n = 1}^\infty  {\frac{{\coth (\pi n) - 1}}{n}} =-0.825861175988396\ldots .
$$
Addendum. Note that, using the Dedekind $\eta$-function,
$$
\sum\limits_{n = 1}^\infty  {\frac{{\coth (\pi n) - 1}}{n}}  =  - 2\log ({\rm e}^{\pi /12} \eta ({\rm i})) = 2\log \left( {\frac{{2\pi ^{3/4} }}{{\Gamma (1/4)}}} \right) - \frac{\pi }{6}
$$
and hence
$$
C = \frac{\pi }{2}\gamma  - G - \frac{{\pi ^2 }}{6} + \pi \log \left( {\frac{{2\pi ^{3/4} }}{{\Gamma (1/4)}}} \right).
$$
