§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)$)
Plot[{z^2 NIntegrate[((4 t)/\[Pi])^(1/2)Exp[-t] EllipticTheta[3, 0, Exp[-t*z^2]]^2, {t, 0, \[Infinity]}],Sum[SquaresR[2, k] z^2/(k*z^2 + 1)^(3/2), {k, 0, 10000}]}, {z, 0, 5}, PlotRange -> {{0, 5}, {0, 25}}]– QCD_IS_GOOD Sep 25 '21 at 12:49