Finding smooth behaviour of infinite sum Define
$$E(z) = \sum_{n,m=-\infty}^\infty \frac{z^2}{((n^2 + m^2)z^2 + 1)^{3/2}} = \sum_{k = 0}^\infty \frac{r_2(k) z^2}{(kz^2 + 1)^{3/2}} \text{ for } z \neq 0$$
$$E(0) = \lim_{z \to 0} E(z) = 2 \pi$$
where $r_2(k)$ is the number of ways of writing $k$ as a sum of two squares of integers. Is $E(z)$ smooth at $z = 0$, and can we evaluate its derivatives $(\partial_z)^n E(z)|_{z = 0}$?

 Footnote 1: There is a physical motivation for these equations, $E(z)$ is the electric field generated by a $2d$ lattice of point charges where the lattice spacing is $z$ and the charge density is held fixed. $E(0)$ here is the continuum limit.  
 Footnote 2: Here is a plot of the behavior of the sum in the second representation I wrote above from $k = 0$ to $k = N$, as a function of $z \in [-2,2]$. The black line is $2 \pi$. It looks like the limit will be very flat at $z = 0$ (possibly all derivatives vanish? That would be cool)



Here's another way to proceed. The Poisson summation formula tells us that summing the function over the lattice is the same as summing its fourier transform. Luckily, we can actually do the fourier transform:
$$f(x) = \int_{\mathbb{R}^2}  \frac{e^{-2 \pi i \langle y , x \rangle}z^2}{(y^2 z^2 + 1)^{3/2}} \mathrm{d}y = \int_\mathbb{R^2} \frac{e^{-2 \pi i v_1 x}z^2}{(v^2 z^2 + 1)^{3/2}} = 2\pi e^{-2\pi |x|/z}$$
where we have used the transformation of coordnates $v_1 = \frac{x_1 y_1 + x_2 y_2}{x}, v_2 = \frac{x_2 y_1 - x_1 y_2}{x}$ to do the first transformation, and then we can integrate over $v_2$ by trigonometric substitution, and over $v_1$ by picking out the residue in the top/bottom half plane of the countour integral.
Using this, we can explicitly write down the function as:
$$E(z) = \sum_{k=0}^\infty  2\pi \ r_2(k) \ e^{-2\pi \sqrt{k}/z} = 2\pi + 8\pi e^{-2\pi/z} + 8\pi e^{-2\sqrt{2}\pi/z} + \dots$$
The 'problem' at $z=0$ in position space has been fixed, and it's clear to see now how all derivatives vanish at $z=0$.
