Closed form for $\sum\limits_{r=-\infty}^\infty \sum\limits_{s=-\infty}^\infty \frac{1}{(k^2\tau^2+(x+2\pi r)^2+(y+2\pi s)^2)^{3/2}}$ I have a problem while trying to find a closed form for the following double sum. 
$$\sum_{r=-\infty}^\infty \sum_{s=-\infty}^\infty \frac{1}{(k^2\tau^2+(x+2\pi r)^2+(y+2\pi s)^2)^{3/2}}$$
I used Mathematica for the evaluation but it just returns the double sum. I am wondering how this kind of infinite sum can be handled.
Any hints or suggestions are welcomed, perhaps by proposing some ideas to simplify it or making some approximations.

@Dr. Wolfgang Hintze Thank you very much for the detailed explanation. I have a question, when I try to plot the double sum I get an error ''NSum :Summand (or its derivative)... - is not numerical at point r = 0 ''? To plot it I just wrote this :

What is the problem in this formula?
 A: This is not a solution, as I haven't found a closed form expression (and I doubt meanwhile that it exists) but a collection of illustrative remarks that might be useful. 
Notice that I have reinserted and simplified the derivation of the Fourier form of the double sum in 6).
1) While the sum
$$f=\sum _{r=-\infty }^{\infty } \sum _{s=-\infty }^{\infty } \frac{1}{\left(k^2 \tau ^2+(2 \pi  r+x)^2+(2 \pi  s+y)^2\right)^{3/2}}$$
is returned unevaluated by Mathematica the corresponding integral has a very simple closed form solution (independent of $x$ and $y$)
$$fi=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\frac{1}{\left(k^2 \tau ^2+(2 \pi  r+x)^2+(2 \pi  s+y)^2\right)^{3/2}}dsdr = \frac{1}{2 \pi  \sqrt{k^2 \tau ^2}}$$
This is of course a purely technical remark as the integral lacks the essential dependence on $x$ and $y$.
2) Numerical investigation shows that but a few summand are sufficient to obtain a reasonable approximation to the complete sum.
3) Plots of the double sum 
3a) 3D-plot for two values of $z = k \tau$
$z = 1$

$z = \pi /2$ (this is the value to be assigned to z acording to the OP)

We see the expected regular spike structure. Where the spikes become steeper with smaller $z$ (notice that the vertical scale is the same in the two pictures).
3b) Cross sections $y = 0$ for different values of z

3c) Maximum value of $f$ as a function of $z$

Notice that $f \simeq 1/z^3$ for small $z$. Which can also be derived directly form the double sum as is shown in section 5).
3d) Minimum value of $f$ as a function of $z$ 

4) Physical interpretation of the double sum.
Consider an infinite array of unit masses placed in the x-y-plane at the points $(x,y) = 2 \pi (r,s)$.
The double sum $f$ multiplied by $z$ is then equal to the vertical component of the gravitational force on a unit mass at the point $(x,y,z)$.
It is clear from this analogue model that the force is strong for small $z$ and decreases for increasing $z$.
5) Asymptotic behaviour 
5a) For small $z$ we have $z \simeq 1/z^3$.
Indeed, for small $z$ we can write
$$f(0,0,z\to 0)\simeq \frac{4}{(2 \pi )^{3/2}}\sum _{r=1}^{\infty } \sum _{s=1}^{\infty } \frac{1}{\left(r^2+s^2\right)^{3/2}}+\frac{1}{z^3}$$
The double sum is a finite constant as can be shown thus
$$fr=\sum _{r=1}^{\infty } \sum _{s=1}^{\infty } \frac{1}{\left(r^2+s^2\right)^{3/2}}<\int _1^{\infty }\int _1^{\infty }\frac{1}{\left(r^2+s^2\right)^{3/2}}dsdr=2-\sqrt{2}\simeq 0.585786\text{...}$$
This proves the statement. 
5b) For $z\to\infty$ we find $f\simeq \frac{1}{2 \pi z}$.    
Indeed, in this case the double sum can be approximated by the double integral as calculated in section 1).
5c) Both limits are easily understood from the physical analogy:     
For small $z$ and $(x,y) = (0,0)$ we are approaching a single mass at the origin. The force we have calculated is equal to the $z \frac{1}{z^3} = \frac{1}{z^2}$, in agreement with Newton's law of gravity.
For large $z$, on the other hand, the grid structure is "invisible" to the point mass, and the force becomes equal to the force of an infinite massive plate, and this is constant. Indeed we have found $f \simeq 1/z$, times $z$ this is constant.
6) Fourier form of the double sum
The double sum $f$, renamed $fc$ here, can also be written as a double Fourier series.
$$fc(x,y,z)=\frac{1}{2 \pi  z}\sum _{n=-\infty }^{\infty } \sum _{m=-\infty }^{\infty } e^{-z \sqrt{m^2+n^2}} \cos (m y) \cos (n x)$$
The derivation starts from the expression obtained in the solution of Pierpaolo Vivo 
$$fe(x,y,z)=\int_0^{\infty } \frac{e^{-t z^2} \vartheta _3\left(\frac{x}{2},e^{-\frac{1}{4\; t}}\right) \vartheta _3\left(\frac{y}{2},e^{-\frac{1}{4\; t}}\right)}{2 \pi ^{3/2} \sqrt{t}} \, dt$$
Subsituting the expression for the Jacobi theta function (see the real part of formula (4) in http://mathworld.wolfram.com/JacobiThetaFunctions.html)
$$\vartheta _3\left(\frac{u}{2},q\right)=\sum _{n=-\infty }^{\infty } q^{n^2} \cos (n u)$$
and doing the $t$-integral gives
$$\cos (m y) \cos (n x) \int_0^{\infty } \frac{\exp \left(-\frac{m^2+n^2}{4 t}-t z^2\right)}{\left(2 \pi ^{3/2}\right) \sqrt{t}} \, dt=\frac{e^{-z \sqrt{m^2+n^2}} \cos (m y) \cos (n x)}{2 \pi  z}$$
The double sum of which is $fc$. Done.
Remark: it is interesting, that replacing the double sum by a double integral gives (hint: transform to polar coordinates)
$$\frac{1}{2 \pi  z}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }e^{-z \sqrt{m^2+n^2}} \cos (m y) \cos (n x)dmdn=\frac{1}{\left(x^2+y^2+z^2\right)^{3/2}}$$
This expression is less trivial than the one in 1) but the double periodic structure is washed out as well.
7) Approximations
We can approximate the function by taking only a finite number $n$ of terms, defined here as letting the summation indices run from $-n$ to $n$.
It turns out that for $z\gtrsim 1$ the Fourier double sum $fc$ requires less terms than the plain sum $f$, while for $z\lesssim 1$ the situation is reversed.
The following graphs show the typical behaviour of the maximum and minimum values as a function of the number of terms.
The case $z = \pi /2$, this is the values of the OP


The case $z = 1/2$


We conclude that if $z$ takes on the value $\pi /2$ of the OP, then the Fourier double sum truncated at index $n = 4$ gives a satisfactory approximation to $f$. Hence we can state (somewhat boldly) that this is the requested closed form expression.
A: Not really a full answer, but too long for a comment. You can always use the identity
$$
\frac{1}{Z^{3/2}}=\frac{2}{\sqrt{\pi}}\int_0^\infty d\xi\ \xi^{1/2}e^{-\xi Z}\ ,
$$
for positive $Z$. Your double sum then reads
$$
\frac{2}{\sqrt{\pi}}\int_0^\infty d\xi\ \xi^{1/2}e^{-\xi k^2\tau^2}\sum_{r=-\infty}^\infty e^{-\xi (x+2\pi r)^2}\sum_{s=-\infty}^\infty e^{-\xi (y+2\pi s)^2}
$$
and evaluating the two sums in terms of Elliptic theta functions, one gets
$$
\frac{2}{\sqrt{\pi}}\int_0^\infty d\xi\ \xi^{1/2}e^{-\xi k^2\tau^2}\frac{\vartheta _3\left(\frac{x}{2},e^{-\frac{1}{4 \xi}}\right)}{2 \sqrt{\pi } \sqrt{\xi}}\frac{\vartheta _3\left(\frac{y}{2},e^{-\frac{1}{4 \xi}}\right)}{2 \sqrt{\pi } \sqrt{\xi}}
$$
$$
=\int_0^{\infty } \frac{e^{-k^2 \xi  \tau ^2} \vartheta _3\left(\frac{x}{2},e^{-\frac{1}{4 \xi }}\right) \vartheta _3\left(\frac{y}{2},e^{-\frac{1}{4 \xi }}\right)}{2 \pi^{3/2}  \sqrt{\xi }} \, d\xi.
$$
The integral does not seem to be easy to compute in closed form, but might be better suited for an asymptotic/approximate evaluation. 
A: @Dr. Wolfgang Hintze, Thank you very much for the contribution. 
The point of this answer is just to propose something which could complete the derivation of $f_{s}$ you proposed in the previous answer. At this step I don't think there is a closed form solution, but a reasonable approximation can be made.
Intuitively, $|\cos(nx) \cos(y n)|$ is spread out equally between 0 and 1. If we approximate it as a constant $c>0$, then the sum become
$f_{s} \simeq \frac{2c}{\pi z}\sum_{n=1}^\infty \sum_{m=1}^\infty e^{-z\sqrt{m^2+n^2}}$.
Again this expression is not tractable with Mathematica, so I thought that approximation can be made using the Euler–Maclaurin formula.
I start by doing a simple modification on $f_{s}$ :
$f_{s} \simeq \frac{2c}{\pi z}[\sum_{n=0}^\infty \sum_{m=0}^\infty e^{-z\sqrt{m^2+n^2}}-2]$.
In one variable, Euler_Maclaurin formula state that $\sum_{n=0}^\infty f(n)=\frac{1}{2}f(0)+\int_{0}^{\infty}f(n) dn$, I will use this twice in what follow.
We know that
$\frac{\pi}{2z^2} = \int_{0}^{\infty}\int_{0}^{\infty}e^{-z\sqrt{m^2+n^2}} dn dm$
then
$\frac{\pi}{2z^2} = \int_{0}^{\infty}[\sum_{m=0}^\infty e^{-z\sqrt{m^2+n^2}}-\frac{1}{2}e^{-zm}]dm$
$\frac{\pi}{2z^2} = \sum_{m=0}^\infty \sum_{n=0}^\infty e^{-z\sqrt{m^2+n^2}}-\frac{1}{2}\sum_{n=0}^\infty e^{-zn}+\frac{1}{4}$
$\frac{\pi}{2z^2} = \sum_{m=0}^\infty \sum_{n=0}^\infty e^{-z\sqrt{m^2+n^2}}-\frac{1}{2} \frac{e^{z}}{e^{z}-1}+\frac{1}{4}$
So, we can now conclude that 
$\sum_{m=0}^\infty \sum_{n=0}^\infty e^{-z\sqrt{m^2+n^2}}= \frac{\pi}{2z^2}+\frac{1}{2} \frac{e^{z}}{e^{z}-1}-\frac{1}{4}.$
This all what I was able to do, but I'am not sure that my assumption ($|\cos(nx) \cos(y n)|$=c) is reasonable to guarantee a good approximation. What do you think ?
