Double summation

Question

Can Mathematica compute a double sums such as:

Sum[If[m == n == 0, 0, 
  1/(m^2 + n^2)^3], {m, -∞, ∞}, {n, -∞, \
∞}]

The analytic expression can be found in https://mathworld.wolfram.com/DoubleSeries.html in terms of Zeta functions, but I am curious why Mathematica does not handle this directly.

In case you are curious, the above sum evaluates to:

1/8 \[Pi]^3 Zeta[3]

Answer 1

As a starting point, I suggest to look at this post

https://math.stackexchange.com/questions/1108246/double-sum-and-zeta-function

We can at least verify numerically that the sum is equal to

$$S(3)=\frac1{\Gamma(3)}\int_0^{\infty}\!\! t^2 \big(\theta_3(0,e^{-t})^2-1\big)\, \mathrm{d}t$$

With $\Gamma(3)=2$, we have

NIntegrate[t^2(EllipticTheta[3,0,Exp[-t]]^2-1),{t,0,Infinity}]/2

NSum[If[m == n == 0, 0, 1/(m^2 + n^2)^3], 
      {m, -Infinity, Infinity}, {n, -Infinity, Infinity}]

(*4.65891*)

Note, it is not easy to quickly evaluate the sum numerically, but with the integral representation it evaluates instantly.

Unfortunately MA does not know the general coefficient of the Lambert series expansion of $\theta_3(0,q)^2$. Thus, it cannot assist in derivation along the lines of the linked math.SE post.

However, if we do manual Lambert series expansion of $$ \theta_3(0,q)^2=1+4\sum_{n=1}^\infty \frac{q^n}{1+q^{2n}}, $$ and subsequently expand $$ \theta_3(0,q)^2-1=4\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m q^nq^{2mn}, $$ MA is able to do the remaining Melin transform and a double sum.

4/Gamma[3] Sum[(-1)^m Integrate[t^2 q^n q^(2m n)/.q->Exp[-t],
                       {t,0,Infinity}],
                        {n,1,Infinity},{m,0,Infinity}]

(* Zeta[3](Zeta[3,1/4]-Zeta[3,3/4])/16 *)

We may want to verify the prefactor of the zeta-function

FullSimplify[(Zeta[3,1/4]-Zeta[3,3/4])/16]

(* Pi^3/8 *)

General solution

If we can establish that $$S(s)\equiv\sum_{i\neq j}\frac{1}{(i^2+j^2)^s}= \frac1{\Gamma(s)}\int_0^{\infty}\!\! t^{s-1} \big(\theta_3(0,e^{-t})^2-1\big)\, \mathrm{d}t\\ =\frac{4}{\Gamma(s)}\sum_{n=1}^{\infty} \sum_{m=0}^{\infty} \int_0^{\infty}\!\! t^{s-1} \exp\big\{-n(2m+1)t\big\}\,\mathrm{d}t $$ MA quickly returns the answer

4/Gamma[s]Sum[(-1)^m Integrate[t^(s-1) q^n q^(2m n)/.q->Exp[-t],{t,0,Infinity}],{n,1,Infinity},{m,0,Infinity}]