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