# Closed form sum for the following series on the euclidean grid.

I am trying to find a closed form solution for the following series. The $$\sqrt{i^2 + j^2}$$ in the exponent comes from distances on the euclidean grid from the origin.

$$x = \sum_{i,j} e^{-\sqrt{i^2 + j^2}}$$

where $$i,j$$ range from $$0$$ to infinity.

It appears this expression is not a geometric series, so I have trouble analyzing it. I did some quick simulations to realize that the value converges quickly. For $$i,j$$ in range (0,40), and using double-precision floating point, the value converges to $$2.95878712840391$$. Altering the range of $$i,j$$ no longer changes the sum because the incremental values are beyond the precision of the floating point decimal.

I would greatly appreciate some help in approaching this series, and if there is a way to represent it in a closed form. Or if there is a way to approximate the answer to a desired precision.

By letting $$r_2(n)=\left|\{(a,b)\in\mathbb{Z}^2:a^2+b^2=n\}\right|$$ we have $$r_2(n) = 4(\chi_4 * 1)(n) = 4\sum_{d\mid n}\chi_4(d)$$, with $$\chi_4(n)=1$$ if $$n\equiv 1\pmod{4}$$, $$\chi_4(n)=-1$$ if $$n\equiv -1\pmod{4}$$ and $$\chi_4(n)=0$$ if $$n$$ is even (i.e. $$\chi_4$$ is the non-principal Dirichlet character $$\!\!\!\pmod{4}$$).
It follows that $$\sum_{i,j\geq 0}e^{-\sqrt{i^2+j^2}} = 1+4\sum_{n\geq 1}(\chi_4*1)(n) e^{-\sqrt{n}}$$ where $$(\chi_4*1)(n)$$ has a moderately erratic behaviour, but $$e^{-\sqrt{n}}$$ converges to zero really fast, such that in order to compute the LHS up to $$N$$ figures it is enough to compute $$\sum_{n=1}^{N^2}(\chi_4*1)(n) e^{-\sqrt{n}}$$. As an alternative, we have $$e^{-\sqrt{n}}=\frac{1}{\sqrt{\pi}}\int_{0}^{+\infty}e^{-s^2/4}e^{-n/s^2}\,ds=\frac{1}{2\sqrt{\pi}}\int_{0}^{+\infty}e^{-\frac{1}{4s}}e^{-ns}\,\frac{ds}{s^{3/2}}$$ hence by setting $$\Theta(x)=\sum_{n\geq 0} x^{n^2}$$ we have $$\sum_{i,j\geq 0}e^{-\sqrt{i^2+j^2}}=\frac{1}{\sqrt{\pi}}\int_{0}^{+\infty}\Theta^2(e^{-s^2}) e^{-\frac{1}{4s^2}}\frac{ds}{s^{2}}$$ where the RHS can be approximated through numerical integration algorithms, especially if combined with the functional identity for the Jacobi $$\Theta$$ function, which is a consequence of the Poisson summation formula. My computations point towards an approximated value of $$\color{green}{2.9587871284039}\color{red}{3}$$.