A computationally attractive integral representation of (even more general) $$\Phi_k(\alpha)=\sum_{n_1,\ldots,n_k\in\mathbb{Z}}e^{-\alpha\sqrt{n_1^2+\ldots+n_k^2}}\qquad(\alpha>0)$$ may be obtained using the following integral (a Laplace transform): $$e^{-\sqrt{s}}=\frac1{2\sqrt\pi}\int_0^\infty t^{-3/2}e^{-st-1/(4t)}\,dt.$$
We get, after the substitution $t=\pi/(\alpha x)^2$,
\begin{align*}
\Phi_k(\alpha)
&=\frac1{2\sqrt\pi}\int_0^\infty t^{-3/2}e^{-1/(4t)}\left(\sum_{n\in\mathbb{Z}}e^{-\alpha^2 n^2 t}\right)^k\,dt
\\&=\frac\alpha\pi\int_0^\infty e^{-(\alpha x)^2/(4\pi)}\vartheta(1/x)^k\,dx,
\end{align*}
where $\vartheta(x):=\sum_{n\in\mathbb{Z}}e^{-\pi(nx)^2}$ is a "theta function", with known $\vartheta(1/x)=x\vartheta(x)$.
Here is a numerical experiment using PARI/GP to illustrate the efficiency. In the code
myexp(x)=if(x<-default(realbitprecision),0,exp(x));
mytheta(x)=1+2*suminf(n=1,myexp(-Pi*(n*x)^2));
optheta(x)=if(x<1,mytheta(1/x)/x,mytheta(x));
myint(f,h)=h*(f(0)+suminf(n=1,f(n*h)+f(-n*h)));
term(n)=2^hammingweight(n)*myexp(-sqrt(norml2(n)));
foo(h)=myint(x->myexp(x-myexp(2*x)/4/Pi)*optheta(myexp(-x))^3,h)/Pi;
goo(m)={my(r=0);forvec(n=vector(3,i,[0,m]),r+=term(n));return(r)};
the function goo(m)
computes the partial sum $\sum_{-m\leqslant n_1,n_2,n_3\leqslant m}$ of $\phi=\Phi_3(1)$, while the function foo(h)
computes the integral for $\Phi_3(1)$ obtained via the substitution $x=e^y$ from the "theta" integral above, using the (simple) trapezoidal rule with step $h$; the remaining functions are for some sort of "optimization". The computations using foo(h)
give $$\phi=\underbrace{25.39268269326689073}_{h=0.1}\underbrace{60512162955540324466}_{h=0.05}\underbrace{338930111641124451230483735111+}_{h=0.02}$$ (in fact foo(0.02)
gives $104$ correct digits of $\phi$, and foo(0.01)
gives $211$ correct digits).
The computations using goo(m)
are much slower (especially w.r.t. the time taken): $$\phi=\underbrace{25.39268}_{m=20}\underbrace{269326689073}_{m=50}\underbrace{60512162955540324466338}_{m=100}\underbrace{930111641124451230483735111+}_{m=200}$$ (here goo(200)
gives $84$ correct digits of $\phi$).