I had trouble in finding a closed form solution for the following series, so now I am trying to find a good approximation for it. The $\sqrt{i^2 + j^2}$ in the exponent comes from distances on the euclidean grid from the origin.


$x = \sum\_{i,j} e^{-a\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. Setting $a = 1$, 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.



Approximate the following series on the euclidean grid
You're surely right that there cannot be a "closed form" for such a series;
but it can still be approximated to any desired precision.
The defining sum
$$
x = x(a)
 = \sum\_{i=0}^\infty \sum\_{j=0}^\infty \exp\Bigl(-a \sqrt{i^2+j^2}\Bigr)
$$
requires only $O(N^2/a)$ terms to approximate to within $\exp(-N)$,
so it converges quickly enough unless $a$ is small;
but in that case one can use a Laurent series
$x(a) = \frac\pi{2a^2} + \sum\_{l=0}^\infty c\_l a^l$
that converges for $a < 2\pi$. The "small" and "large" ranges overlap,
so one can check the analysis by computing $x(a)$ both ways
for some values of $a$ in say $[1/2, 2]$ and noting that
they agree to within the computed precision.


Instead of $x(a)$, which is a sum over integer $(i,j)$
in the first quadrant, it is more natural to sum over all integer pairs:
define
$$
F(a) = \sum\_{i=-\infty}^\infty \sum\_{j=-\infty}^\infty
 \exp \Bigl(-a \sqrt{i^2+j^2}\Bigr).
$$
It is easy to go between $x(a)$ and $F(a)$,
because each term in the $x(a)$ sum appears $4$ times in the $F(a)$ sum
(for $(i,j)$, $(i,-j)$, $(-i,j)$, and $(-i,-j)$),
except when $i=0$ or $j=0$ --- and those terms *do* give geometric series,
so we soon find
$$
F(a) = 4 x(a) - \frac4{1 - e^{-a}} + 1.
$$
We shall prove that for $0 < a < 2\pi$ we have
$$
F(a) = \frac{2\pi}{a^2} + \sum\_{k=1}^\infty c\_{2k-1} a^{2k-1},
$$
where the coefficients $c\_1,c\_3,c\_5,\ldots$ are given by
$$
c\_{2k-1} = 4 \left({-3/2 \atop k-1}\right) (2\pi)^{-2k}
 \zeta\bigl(k+\frac12\bigr) L\bigl(k+\frac12, \chi\_4\bigr).
$$
Here $\zeta$ is the Riemann zeta function and 
$L(\cdot, \chi\_4)$ is the Dirichlet $L$-function
$$
L(s, \chi\_4) = \sum\_{n=0}^\infty \frac{(-1)^n}{(2n+1)^s}
 = 1 - \frac1{3^s} + \frac1{5^s} - \frac1{7^s} + - \cdots.
$$
These factors of $c\_{2k-1}$ are not elementary but can be
efficiently approximated to any desired accuracy,
so the same is true of the Laurent series. Here is a **gp** implementation:



```
L4(s) = sumalt(n=0,(-1)^n/(2*n+1)^s)

{
F(a) = 2*Pi/a^2 + 8*Pi * suminf(k=1,
   binomial(-3/2,k-1) * (2*Pi)^(-2*k-1) * zeta(k+1/2) * L4(k+1/2) * a^(2*k-1))
}

```

To obtain the series expansion of $F(a)$, we first apply
the 2-dimensional
[Poisson
 summation formula](https://en.wikipedia.org/wiki/Poisson_summation_formula)
$$
\sum\_{i=-\infty}^\infty \sum\_{j=-\infty}^\infty f(i,j)
= \sum\_{m=-\infty}^\infty \sum\_{n=-\infty}^\infty \hat f(m,n),
$$
where $\,f: {\bf R}^2 \to \bf C$ is "any" function, and
$\,\hat f$ is its Fourier transform
$$
\hat f(y\_1,y\_2) = \mathop{\iint}\_{(x\_1,\,x\_2) \in {\bf R}^2}
 e^{-2\pi i (x\_1 y\_1 + x\_2 y\_2)} f(x\_1,x\_2) \, dx\_1 \, dx\_2.
$$
For $f(x\_1,x\_2) = \exp\bigl( -a \sqrt{x\_1^2+x\_2^2} \bigr)$
the Fourier transform is elementary, though it takes a bit of work
(write the double integral in polar coordinates,
and integrate with respect to $r$ first). We find
$$
F(a) = 2 \pi a \sum\_{m=-\infty}^\infty \sum\_{n=-\infty}^\infty
 \bigl(a^2 + 4\pi^2(m^2+n^2)\bigr)^{-3/2}.
$$
This sum converges, but too slowly to use as it stands.
(The singularity of $\,f$ at $(0,0)$ makes $\,\hat f$ decay slowly
at infinity, though the decay is still fast enough to justify
the use of Poisson summation.) So for each $(m,n) \neq (0,0)$
we write $(a^2 + 4\pi^2(m^2+n^2))^{-3/2}$ as a binomial expansion
$$
(2\pi)^{-3} (m^2+n^2)^{-3/2} \sum\_{k=1}^\infty
 \left({-3/2 \atop k-1}\right) \bigl(a^2 / (4\pi^2 (m^2+n^2))\bigr)^{k-1},
$$
which converges if $|a| < 2 \pi$, and then combine terms with
$a^2$ appearing to the same power. The formula for the coefficients
$c\_{2k-1}$ then follows using the known factorization
$$
\mathop{\sum\!\sum}\_{\scriptstyle(m,n) \in {\bf Z}^2 \atop \scriptstyle(m,n) \neq (0,0)}
 (m^2+n^2)^{-s} = 4 \zeta(s) \, L(s,\chi\_4)
$$
which encodes unique factorization in ${\bf Z}[i]$.


