# Convergence check of a multiple series

I want to check that the following series $$s_1=\sum_{i,j=1}^{\infty}\frac{ij}{(i+j)^{a/2+1}}$$ diverges for $$a\le 6$$, where $$a$$ is a positive integer. I tried with SumConvergence: SumConvergence[(i j)/(i + j)^(a/2 + 1), {i, j}]without getting an answer. Can anyone help?

The simplest approach is to replace the sum by an integral:

Integrate[(x y)/(x + y)^(a/2 + 1), {x, 1, ∞}, {y, 1, ∞}]

(*... Re[a] > 6 *)


For greater clarity, I will consider a series

 Sum[(i j)/(i + j)^k, {i, 1, Infinity}]


Mathematica does not generally sum this series, but we can derive a general formula

 Table[Sum[(i j)/(i + j)^k, {i, 1, Infinity}], {k, 4, 10}]
(* {-(1/6) j (3 PolyGamma[2, 1 + j] + j PolyGamma[3, 1 + j]),
1/24 j (4 PolyGamma[3, 1 + j] + j PolyGamma[4, 1 + j]),
-(1/120) j (5 PolyGamma[4, 1 + j] + j PolyGamma[5, 1 + j]),
1/720 j (6 PolyGamma[5, 1 + j] + j PolyGamma[6, 1 + j]),
-((j (7 PolyGamma[6, 1 + j] + j PolyGamma[7, 1 + j]))/5040),
(j (8 PolyGamma[7, 1 + j] + j PolyGamma[8, 1 + j]))/40320,
-((j (9 PolyGamma[8, 1 + j] + j PolyGamma[9, 1 + j]))/362880)} *)


From here we have the same result

 Table[(-1)^(k + 1) 1/(k - 1)! j ((k - 1)*PolyGamma[k - 2, 1 + j]
+ j PolyGamma[k - 1, 1 + j]), {k, 4, 10}]


Now we will see how these expressions behave for large j

 Table[Series[(-1)^(k + 1) 1/(k - 1)! j ((k - 1)*
PolyGamma[k - 2, 1 + j] + j PolyGamma[k - 1, 1 + j]), {j,
Infinity, k}], {k, 4, 10}]
(* {1/(6 j)-1/(12 j^3)+O((1/j)^5),
1/(12 j^2)-1/(12 j^4)+O((1/j)^6),
1/(20 j^3)-1/(12 j^5)+O((1/j)^7),
1/(30 j^4)-1/(12 j^6)+O((1/j)^8),
1/(42 j^5)-1/(12 j^7)+O((1/j)^9),
1/(56 j^6)-1/(12 j^8)+O((1/j)^10),
1/(72 j^7)-1/(12 j^9)+O((1/j)^11)} *)


The terms with the highest weight are

 Table[1/((k - 1)*(k - 2)*j^(k - 3)), {k, 4, 10}]
(* {1/(6 j), 1/(12 j^2), 1/(20 j^3), 1/(30 j^4), 1/(42 j^5), 1/(56 j^6), 1/(72 j^7)} *)


The sum_{j=1..infinity} 1/j^n converges for n>1, so we have

  k - 3 > 1 /. k -> a/2 + 1 // Simplify
(* a > 6 *)

\$Version

(* "13.3.1 for Mac OS X ARM (64-bit) (July 24, 2023)" *)

Clear["Global*"]

sum[α_] := sum[α] =
Sum[i*j/(i + j)^(α/2 + 1), {i, 1, ∞}, {j, 1, ∞}]

Table[{α, sum[α]}, {α, 5, 10}]


sum is only defined for even values greater than 6

data = Table[{α, sum[α]}, {α, 8, 30, 2}]

(* {{8, -(1/540) π^2 (-15 + π^2)}, {10, 1/6 (Zeta[3] - Zeta[5])}, {12, (
21 π^4 - 2 π^6)/11340}, {14,
1/6 (Zeta[5] - Zeta[7])}, {16, -((π^6 (-10 + π^2))/56700)}, {18,
1/6 (Zeta[7] - Zeta[9])}, {20, (π^8 (99 - 10 π^2))/5613300}, {22,
1/6 (Zeta[9] - Zeta[11])}, {24, (π^10 (6825 - 691 π^2))/
3831077250}, {26,
1/6 (Zeta[11] - Zeta[13])}, {28, (π^12 (691 - 70 π^2))/
3831077250}, {30, 1/6 (Zeta[13] - Zeta[15])}} *)


For even values the sums are

f[α_] = 1/6 (Zeta[(α - 4)/2] - Zeta[α/2]);


Verifying the sums for even values,

And @@ Table[sum[α] == f[α], {α, 8, 30, 2}]

(* True *)


f has a singularities at α == 2 and α == 6

FunctionDomain[f[α], α]

(* α < 2 || 2 < α < 6 || α > 6 *)

Plot[f[α], {α, 0, 10}, PlotRange -> All]


Show[
Plot[f[α], {α, 7, 30},
PlotRange -> All],
DiscretePlot[sum[α], {α, 8, 30, 2},
PlotStyle -> Red]]


Show[
LogPlot[f[α], {α, 7, 30},
PlotRange -> All],
DiscretePlot[sum[α], {α, 8, 30, 2},
PlotStyle -> Red,
ScalingFunctions -> "Log"]]


Use f for the sum, it is much faster

sum[48] // AbsoluteTiming

(* {26.9225, (π^22 (2332820490 - 236364091 π^2))/1211517431782539131250} *)

f[48] // AbsoluteTiming

(* {0.000031,
1/6 ((155366 π^22)/13447856940643125 - (236364091 π^24)/
201919571963756521875)} *)

%[[2]] == %%[[2]] // Simplify

(* True *)
`