Multiple series calculation

Question

Let $n$ be a positive integer. I would like to find a numerical evaluation of the convergent (!) series $$ S_{n,s}=\sum_{k\in \mathbb Z^{n}}\frac{1}{(1+\vert k\vert^{2})^{s/2}},\quad s> n, $$ where $\vert k\vert$ stands for the standard Euclidean norm in $\mathbb R^n$. As mentioned above, the convergence is easily proven and follows from a comparison with an integral. However the numerical evaluation requires to take care of the small values of $\vert k\vert$. Maybe there is a close formula either for $S_{n,s}$ or for $$ \tilde S_{n,s}=\sum_{k\in \mathbb Z^{n}\backslash\{0\}}\vert k\vert^{-s}. $$ I am in fact interested in the case $n\ge 2$.

Answer 1

$\newcommand{\vpi}{\varphi} \newcommand{\vp}{\varepsilon}$ We have \begin{align*} S_{n,s}&=\sum_{k\in \mathbb Z^{n}}\frac{1}{\Gamma(s/2)}\int_0^\infty u^{s/2-1} e^{-(1+|k|^2)u}\,du \\ &=\frac{1}{\Gamma(s/2)}\,\int_0^\infty du\, u^{s/2-1}e^{-u} \sum_{k\in \mathbb Z^{n}}e^{-|k|^2u} \\ &=\frac{1}{\Gamma(s/2)}\,\int_0^\infty du\, u^{s/2-1}e^{-u} \Big(\sum_{k\in \mathbb Z}e^{-k^2u}\Big)^n \\ &=\frac{1}{\Gamma(s/2)}\,\int_0^\infty du\, u^{s/2-1}e^{-u}\, \vpi(u)^n, \end{align*} where $\vpi(u):=\vartheta_3(0,e^{-u})$ and $\vartheta_a(u,q)$ is the elliptic theta function.

So, instead of an $n$-fold sum, we now only need to compute an ordinary integral of a special function.

---

E.g., it takes the Mathematica command NSum[] 22 sec to output the value $12.6028$ for the double sum $S_{2,2.5}$ computed according to its definition. In comparison, it takes Mathematica command NIntegrate[] only 0.054 sec to output the value $12.6001$ for the integral expression of this double sum, and only 0.023 sec to output the value $22.0645$ for the integral expression of the triple sum $S_{3,3.5}$. (I have not attempted an evaluation of a triple sum using NSum[]. Even for $n=2$, I tend to believe NIntegrate[] more than NSum[], which latter lacks in accuracy and speed even for ordinary sums.)

Here is a copy of the corresponding Mathematica notebook:

---

Added: As pointed out in the comments and answer by Henri Cohen, the Mathematica calculations of the above integral expression for $S_{n,s}$ are probably faulty. This can be helped as follows. Note that the series \begin{equation} \vpi(u)=\sum_{k\in \mathbb Z}e^{-k^2u} \end{equation} converges very fast if $u>0$ is bounded away from $0$. In particular, for any natural $m$ and all $u\ge\pi$ \begin{equation} \vpi_m(u)<\vpi(u)<\vpi_m(u)(1+\vp_m), \end{equation} where \begin{equation} \vpi_m(u):=1+2\sum_{k=1}^{m-1}e^{-k^2u},\quad \vp_m:=\frac{2e^{-m^2\pi}}{1-e^{-5\pi}}<2.000001e^{-m^2\pi}; \end{equation} here we used the inequality $\sum_{k=m}^\infty e^{-k^2u}\le\sum_{k=m}^\infty e^{-k^2\pi}$ and bounded the latter sum by the sum of a geometric series.

Using the Jacobi identity $\vpi(u)=\pi^{1/2}u^{-1/2}\vpi(\pi^2/u)$ for $u>0$ (which is an instance of the Poisson summation formula), we now represent $\vpi(u)$ as the sum of a series that converges very fast if $u>0$ is bounded away from $\infty$. In particular, for any natural $m$ and all $u\in(0,\pi]$ \begin{equation} \pi^{1/2}u^{-1/2}\vpi_m(\pi^2/u)<\vpi(u)<\pi^{1/2}u^{-1/2}\vpi_m(\pi^2/u)(1+\vp_m). \end{equation}

Thus, \begin{equation} S_{m;n,s}<S_{n,s}<S_{m;n,s}(1+\vp_m), \tag{$\ast$} \end{equation} where \begin{align*} S_{m;n,s}&:=\frac{1}{\Gamma(s/2)}\,\int_\pi^\infty du\, u^{s/2-1}e^{-u}\vpi_m(u)^n \\ &+\frac{\pi^{n/2}}{\Gamma(s/2)}\,\int_0^\pi du\, u^{s/2-1-n/2}e^{-u}\vpi_m(\pi^2/u)^n, \end{align*} so that now we only have to integrate elementary functions.

E.g., using $(\ast)$ with $m=2$, in 0.014 sec we get the approximate value $12.602759835$ for $S_{2,2.5}$, which agrees with Henri Cohen's Pari/GP calculation.

Answer 2

What follows is an answer to Iosif Pinelis, hence not directly to the OP, but may be useful nonetheless.

1) Although the elliptic theta function is available in GP under the name mfTheta(), I just discovered a bug (to be corrected this week-end), so I will not use it.

2) The correct way to sum in a stupid way (i.e., without using specific clever summation methods) in GP is the suminf command which stops when the summand becomes negligible, hence is useable only for geometrically convergent series, not for $1/n^2$. Hence one could in principle write th(u)=1+2*suminf(n=1,exp(-n^2*u)) but this would fail for two different reasons. First, the convergence for $u$ close to $0$ is too slow: for this, we use the functional equation. Second, when $u$ is large, even exp(-u) will underflow. Thus I first write the following:

th0(u) = if(u>1000, 1, 1+2*suminf(n=1, exp(-n^2*u)));

th(u) = if(u>Pi, th0(u), sqrt(Pi/u)*th0(Pi^2/u));

3) The integration process: the most robust method is the doubly-exponential integration method, which in GP is the command intnum. But we must be careful about the endpoints:

-- At infinity, the integrand is like exp(-u) times something much less important, so is coded [oo,1] (exponential decrease).

-- At 0, the integrand is asymptotic to a constant times u^{s/2-1-n/2}, the n/2 comming from theta.

Thus, the GP command that I use is:

S(n,s) = 1/gamma(s/2)*intnum(u=[0,s/2-1-n/2], [oo,1], u^(s/2-1)*exp(-u)*th(u)^n);

4) Concerning accuracy of integration: you are right, nothing is guaranteed. However, consider for instance S(2,2.5). Let $S_1$, $S_2$, $S_3$, and $S_4$ be the result at $38$, $77$, $154$, and $308$ decimals respectively (note that even at $154$ digits the computation requires only $0.3$ seconds). Then $S_1-S_4$, $S_2-S_4$, $S_3-S_4$ are of the order of $10^{-25}$, $10^{-37}$, $10^{-59}$, admittedly not great, but with extreme confidence I would say that $S_1$ is correct to $24$ decimal digits, $S_2$ to $36$, and $S_3$ to $58$, thus giving the answers that I wrote in my comments.