This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 1

The virial series of the dD ideal Bose gas p(ρ)=(∑_(k=1)^∞)(Subscript[c, k]^((d)))(ρ^k) is defined implicitly through the two series p(z)=(∑_(k=1)^∞)(Subscript[α, k]^((d)))(z^k) and ρ(z)=(∑_(k=1)^∞)k (Subscript[α, k]^((d)))(z^k) where (Subscript[α, k]^((d)))=(k^(-d/2-1))(here p is the pressure, ρ the density, and z the fugacity) [✶], [✶], [✶], [✶], [✶]. For d=4, the series p(ρ) is conjectured to have the convergence radius (π^2)/6 (=ζ(d=4)) [✶]. Calculate a truncated form of p(ρ) and its associated diagonal Padé approximation. The smallest root of the denominator of the Padé rational function approximates the convergence radius [✶], [✶]. Confirm to three digits that the convergence radius is (π^2)/6.

makepOfρSeries[o_, d_, prec_] := 
Module[{b, zOfρSeries},
 b[k_] := k^(-d/2 - 1);
 zOfρSeries = InverseSeries[N[Sum[k b[k] z^k, {k, o}], prec] + 
                            O[z]^(o + 1), ρ];
 Normal[Sum[b[k] zOfρSeries^k, {k, o}]]];

While we can in principle get an exact result for the first terms, they grow quickly in size and no obvious pattern for a closed-form of the coefficients comes to mind.

Factor /@ makepOfρSeries[6, d, Infinity]

For even dimensions, the coefficients become rational numbers.

makepOfρSeries[12, 4, Infinity]

For explicit numeric coefficients, the series coefficients do not grow in size and we can easily calculate, say, 250 terms of the series for p(ρ).

pOfρSeries = makepOfρSeries[250, 4, 200];
(* do last terms have verified digits? *) Precision[pOfρSeries]

Applying the root test to the so-calculated series coefficients shows that the convergence radius is a number near (π^2)/6≈1.64493….

kRoots = MapIndexed[If[#1 != 0, {#2[[1]], Abs[#1]^(-1/#2[[1]])}, 
                     Sequence @@ {}]&, CoefficientList[pOfρSeries, ρ]];
ListPlot[kRoots]

SequenceLimit[Drop[Last /@ kRoots, 50], WynnDegree -> 10]

Assuming a singularity of the form (((ρ-(SuperStar[ρ])))^α), we can use the Domb–Sykes interpolation [✶] to approximate the singularity (SuperStar[ρ]). Taking the last 20 series terms into account, this results again in a value for (SuperStar[ρ]) that is different from (π^2)/6 by only 0.02%.

clρ = CoefficientList[pOfρSeries, ρ];
pairs = Partition[Take[MapIndexed[{#2[[1]], #1}&, 
                                  Rest[clρ]], All], 2, 1];
ratios = {1/#[[2, 1]], #[[2, 2]]/#[[1, 2]]}& /@ pairs;
(* interpolate term[n]/term[n - 1] as a function of 1/n *)
1/Coefficient[Fit[ratios, {1, x}, x], x, 0] - Pi^2/6 // N

Now, we form the [100/100] diagonal Padé approximation from the series pOfρSeries and extract its denominator.

Needs["Calculus`Pade`"]

pOfρApproxDen[ρ_] = Pade[pOfρSeries /. c_ ρ^k_ :> SetPrecision[c, 400] ρ^k, 
                         {ρ, 0, 100, 100}] // Denominator;

Calculating the smallest root gives an approximation for the convergence radius that deviates from (π^2)/6 by only 0.04%.

roots = Roots[SetPrecision[N[pOfρApproxDen[ρ], 120], Infinity] == 0, ρ];

Min[Abs[Cases[N[roots, 40], _?NumberQ, Infinity]]] - Pi^2/6