My question is concerned with the convergence of the sum
$$\displaystyle \sum_{\substack{\ \ m \in \mathbb{Z}^d} \\ {\ \ \ m \neq 0}} \frac{1}{\|m\|_{d}^s},$$
where $d \in \mathbb{N}$, and $\|m\|_d := \sqrt{|m_1|^2 + ... + |m_d|^2}$ denotes the usual Euclidean norm on $\mathbb{R}^d$. Is it true that this sum converges whenever $\mathrm{Re}(s) > d$?
Obviously, when $d = 1$, we get the usual definition of the Riemann zeta function, which we know converges for $\mathrm{Re}(s) > 1$ and admits an analytic continuation to the rest of the complex plane. The convergence of similar sums has been investigated -- for instance, see this question on the convergence of the Riemann zeta function in higher dimensions, and there are similar results out there for modified zeta functions such as the so-called Barnes zeta function. These two results intuitively seem to suggest that one might expect the sum to converge for $\mathrm{Re}(s) > d$. Mathematica does not seem to be able to deal with the sum for $d \geqslant 2$.
Can anyone confirm if this is true (or if it is false, provide conditions under which this sum does converge)? Moreover, does anyone know of any reference where this question has been discussed?
Yes your series converges for $Re(s) > d$ and admits a meromorphic continuation to the whole complex plane.
Let $$\zeta_d(s) = \sum_{ m \in (\mathbb{Z}^d)^*} \|m\|^{-s} = \sum_{ m \in (\mathbb{Z}^d)^*} (\sum_{k=1}^d m_k^2)^{-s/2}$$
Then using as usual $\pi^{-s/2}\Gamma(s/2) a^{-s} = \int_0^\infty x^{s/2-1} e^{-\pi a^2 x}dx$ (change of variable $y = \pi a^2 x$)
you get $$\Lambda(s) = \pi^{-s/2}\Gamma(s/2)\zeta_d(s) = \int_0^\infty x^{s/2-1}\sum_{ m \in (\mathbb{Z}^d)^*} e^{-\pi\|m\|^2x}dx$$
Now letting $\theta_d(x) = 1+\sum_{ m \in (\mathbb{Z}^d)^*} e^{-\pi\|m\|^2x} = \sum_{ m \in \mathbb{Z}^d} e^{-\pi\|m\|^2x}$, note that $f(m) = e^{-\pi \|m\|^2 } =\prod_{k=1}^d e^{-\pi m_k^2 }$ is its own $d$-dimensional Fourier transform, so the Poisson summation formula applies to give $$\theta_d(1/x) = \sum_{ m \in \mathbb{Z}^d} f(\frac{m}{\sqrt{x}}) = \sum_{ n \in \mathbb{Z}^d} x^{d/2}\hat{f}(n\sqrt{x}) = x^{d/2} \theta(x)$$ and hence $$\int_0^1 x^{s/2-1}(\theta_d(x)-1)dx \overset{y = 1/x}= \int_1^\infty y^{-s/2-1}(\theta_d(1/y)-1)dy$$ $$ = \int_1^\infty y^{-s/2-1}(y^{d/2}\theta_d(y)-1)dy = \frac{2}{s-d}-\frac{2}{s}+\int_1^\infty y^{-s/2-1}y^{d/2}(\theta_d(y)-1)dy $$ i.e. $$\Lambda(s)= \int_0^1+\int_1^\infty x^{s/2-1} (\theta_d(x)-1)dx = \frac{2}{s-d}-\frac{2}{s}+\underbrace{\int_1^\infty (x^{s/2-1}+x^{(d-s)/2-1}) (\theta_d(x)-1)dx}_{\text{entire since } \theta_d(x)-1 \, = \,\mathcal{O}(e^{-x})} $$
That is meromorphic on the whole complex plane with a pole of order $1$ at $s=d$ and $s=0$, and $\Lambda(s) = \Lambda(d-s)$.
Q.e.d the pole at $s= d$ shows your series converges for $Re(s) > d$, and that it admits an analytic continuation and a functional equation.