# Double summation

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]

• Using Piecewise instead of If doesn't seem to work either: Piecewise[{{0, n == m == 0}, {1/(m^2 + n^2)^3, True}}] – flinty Jun 24 '20 at 20:38
• This post might be interesting Double series over primes – Artes Jun 24 '20 at 23:41
• I don't find it even a weakness. The statement "The analytic expression can be found in mathworld.wolfram.com/DoubleSeries.html in terms of Zeta functions" does not correspond to reality. Up to mathworld.wolfram.com/DoubleSeries.html , the sum under consideration is expressed as as the sum of a certain series times by the sum of another series. – user64494 Jun 25 '20 at 6:55
• I am not sure if I follow what you said, user64494. If you look at expression (38) in mathworld.wolfram.com/DoubleSeries.html tells you that the above sum is simply 1/8 [Pi]^3 Zeta[3]''' – user12588 Jun 25 '20 at 13:26
• Yes, iif $s=3$ , then $\beta(3)=\frac {\pi^3} 8$.up to mathworld.wolfram.com/DirichletBetaFunction.html . I was not wright in the above comment and I am sorry for that. However, I still don't find it even a weakness. I think this is art for art's sake. NSum produces the correct numeric result. – user64494 Jun 25 '20 at 14:25

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}]
`