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)



§1. Basic results
We have to investigate the sum
$$E_s(z) = \sum_{n,m=-\infty}^\infty \frac{z^2}{((n^2 + m^2)z^2 + 1)^{3/2}}\tag{1.1}$$
Writing
$$\frac{z^2}{((n^2 + m^2)z^2 + 1)^{3/2}}=\int_{0}^{\infty} z^2 \sqrt{\frac{4 t}{\pi} }e^{-t(z^2 (m^2 + n^2) + 1)}\,dt\tag{1.2}$$
and doing the sums under the integral the double sum becomes this integral
$$E_i(z) = z^2 \int_0^{\infty}\sqrt{\frac{4 t}{\pi} }e^{-t} \vartheta _3\left(0,e^{-t z^2}\right)^2 \,dt\tag{1.3}$$
where $\vartheta _3$ is a Jacobi theta function.
This expression is better suitable to be plotted (using the NIntegrate command of Mathematica) than the sum.
Here is a picture of your function $E(z)$ for $z>0$ (for $z<0$ we have $E(-z)=E(z)$).

Notice that it is rather different from the plot in the OP.
Further results are:
For $z \to 0$ we have $E(z) = 2 \pi$, for $z>>1$ we have $E(z) \simeq z^2$.
All derivatives of $E(z)$ with respect to $z$ vanish for $z\to 0$ like $z^p e^{-1/z^2}$ (with some power $p$).
§2. Asymptotic expansions
Expanding the integrand $i(z,t)$ with respect to $z$ we obtain

*

*for $z\simeq 0$
$$\begin{align}i = \frac{2\sqrt{\pi }}{\sqrt{t}} e^{-\frac{8 \pi ^2}{t z^2}-t} \left(2 e^{\frac{3 \pi ^2}{t z^2}}+e^{\frac{4 \pi ^2}{t z^2}}+2\right) \times\left(e^{\frac{3 \pi ^2}{t z^2}} \left(2-t z^2\right)+e^{\frac{4 \pi ^2}{t z^2}}+2\right)\end{align}\tag{2.1}$$
integrating over $t$ results in
$$\begin{align}E(0<z\simeq 0)=\pi  \left(-e^{-\frac{2 \pi }{z}} \left(z^2+2 \pi  z-8\right)-2 e^{-\frac{2 \sqrt{2} \pi }{z}} \\
\left(z^2+2 \sqrt{2} \pi  z-4\right)-2 e^{-\frac{2 \sqrt{5} \pi }{z}} \left(z^2+2 \sqrt{5} \pi  z-8\right)\\
+8 e^{-\frac{4 \pi }{z}}+8 e^{-\frac{4 \sqrt{2} \pi }{z}}+2\right)\end{align}\tag{2.2}$$
the leading terms of which are
$$\begin{align}E(0<z\simeq 0)=\pi  \left(2 + 8 e^{-\frac{2 \pi }z}+ 8 e^{-\frac{4 \pi }z}\right)\end{align}\tag{2.3}$$


*for $z\to \infty$
Developing first $\vartheta _3\left(0,e^{-u}\right)^2$ for $u\to \infty$ and then replacing $u \to t z^2$ we find for the integrand
$$i(z\to +\infty) \simeq \frac{2 e^{-t} \sqrt{t} z^2 \left(2 \left(e^{-4 t z^2}+e^{-t z^2}\right)+1\right)^2}{\sqrt{\pi }}\tag{2.4}$$
After integration over $t$ this leads to the interesting result
$$\begin{align}E(z\to \infty) \simeq  z^2 \left(1+\frac{4}{\left(1+z^2\right)^{3/2}}+\frac{4}{\left(1+2 z^2\right)^{3/2}}\\
+\frac{4}{\left(1+4 z^2\right)^{3/2}}+\frac{8}{\left(1+5 z^2\right)^{3/2}}+\frac{4}{\left(1+8 z^2\right)^{3/2}}\right) \end{align}\tag{2.5}$$
Here the following two sequence of numbers appear: $\{0,1,2,4,5,8,...\}$ as a factor before $z^2$ in the denominator and $\{1,4,4,4,8,4,...\}$ in the numerator.
The first series (if prolonged by taking into account higher order terms) is easily identified in OEIS  as
A001481 Numbers that are the sum of 2 squares.
The second series is A004018 Theta series of square lattice (or number of ways of writing n as a sum of 2 squares). Often denoted by r(n) or r_2(n).
Hence the asymptotic form $(2.5)$ leads automatically to the second sum in the OP (the one containing $r_2(k)$)
