(* Manuel Kauers's Asymptotics package. Download from http://www.risc.jku.at/research/combinat/software/Asymptotics/ *)
<<Asymptotics.m

(* Christoph Koutschan's HolonomicFunctions package. Download from http://www.risc.jku.at/research/combinat/software/HolonomicFunctions/ *)
<<HolonomicFunctions.m

HolonomicFunctions package by Christoph Koutschan, RISC-Linz, Version 1.5.1 (09.08.2011)
⟶ Type  ?HolonomicFunctions  for help

Some useful procedures

(* Given a differential equation for f(z), compute the recurrence for the partial sums of the Taylor coefficients of f(z). *)
Clear[RecurrenceForReturnProbability];
RecurrenceForReturnProbability[ann_,Der[z_],f_[n_]]:=
Module[{ann1,ann2,rec},
(* View f(z) as a bivariate function f(n,z) *)
ann1=ToOrePolynomial[Append[ann,S[n]-1],OreAlgebra[Der[z],S[n]]];
ann2=Annihilator[1/(1-z)/z^(n+1),{Der[z],S[n]}];
(* This computes a recurrence for the Taylor coefficients of f(z)/(1-z) via Cauchy's integral formula. *)
rec=FindCreativeTelescoping[DFiniteTimes[ann1,ann2],Der[z]][[1,1]];
Return[ApplyOreOperator[Factor[rec],f[n]]];
];

(* given a recurrence rec in f[n], compute the values {f[0],f[1],...,f[bound]} 
where inits are the initial values {f[0],...,f[d-1]} with d being the order of the recurrence *)
Clear[UnrollRecurrence];
UnrollRecurrence[rec1_,f_[n_],inits_,bound_]:=
Module[{i,x,vals=inits,rec=rec1},
If[Head[rec]=!=Equal,rec=(rec==0)];
rec=rec/.n->n-Max[Cases[rec,f[n+a_.]:>a,Infinity]];
Do[
AppendTo[vals,Solve[rec/.n->i/.f[i]->x/.f[a_]:>vals[[a+1]],x][[1,1,2]]];
,{i,Length[inits],bound}];
Return[vals];
];

(* Given the first values {f[0],...,f[m]} of a sequence f[n] and a basis of its asymptotic solutions,
compute the limit Limit[f[n], n->Infinity]. *)
Clear[MyLimit];
MyLimit[data_List,asym_,n_]:=
Module[{c,d=Length[asym],pos,ansatz,sol},
pos=Length[data]+Range[-d,-1];
ansatz=Array[c,d].asym;
sol=Solve[((ansatz/.n->#)==data[[#+1]])&/@pos,Array[c,d]][[1]];
Return[N[c[d]/.sol,200]];
];

Example: The two-dimensional lattice

ann=Annihilator[1/(1-z*x1*x2)/Sqrt[1-x1^2]/Sqrt[1-x2^2],{Der[x1],Der[x2],Der[z]}]

FindCreativeTelescoping[ann,{Der[x1],Der[x2]}]

OreReduce[%[[1,1]]+Der[x1]**%[[2,1,1]]+Der[x2]**%[[2,1,2]],Together[ann],Extended->True]

FindCreativeTelescoping[ann,Der[x1]]

FindCreativeTelescoping[First[%],Der[x2]]

ApplyOreOperator[%[[1,1]],P[z]]==0

DSolve[{%,P[0]==1,P'[0]==0},P[z],z]

The four-dimensional face-centred cubic lattice

Differential equation for the LGF

The lattice Green's function is given by the following four-fold integral:

TraditionalForm[HoldForm[Integrate[1/(1-z/6*(Cos[k1]*Cos[k2]+Cos[k1]*Cos[k3]+Cos[k2]*Cos[k3]+Cos[k1]*Cos[k4]+Cos[k2]*Cos[k4]+Cos[k3]*Cos[k4])),{k1,0,Pi},{k2,0,Pi},{k3,0,Pi},{k4,0,Pi}]]]

After the substitutions xi=cos(ki) the integrand transforms to:

integrand=1/(1-z/6*(x1*x2+x1*x3+x1*x4+x2*x3+x2*x4+x3*x4))/(Sqrt[1-x1^2] Sqrt[1-x2^2] Sqrt[1-x3^2] Sqrt[1-x4^2]);

ann0=Annihilator[integrand,{Der[x1],Der[x2],Der[x3],Der[x4],Der[z]}]

ann0 is a system of PDEs for the integrand.

Timing[{ann1,delta1}=CreativeTelescoping[ann0,Der[x1],{Der[x2],Der[x3],Der[x4],Der[z]}];]

(* Test *)
Timing[OreReduce[MapThread[(#1+Der[x1]**#2)&,{ann1,delta1}],ann0]]

Hence, ann1 is a system of PDEs for the first integral (integration w.r.t. x1).

Timing[{ann2,delta2}=CreativeTelescoping[ann1,Der[x2],{Der[x3],Der[x4],Der[z]}];]

(* Test *)
Timing[OreReduce[MapThread[(#1+Der[x2]**#2)&,{ann2,delta2}],ann1]]

Hence, ann2 is a system of PDEs for the second integral (integration w.r.t. x1 and x2).

Timing[{ann3,delta3}=FindCreativeTelescoping[ann2,Der[x3]];]

(* Test *)
Timing[OreReduce[MapThread[(#1+Der[x3]**First[#2])&,{ann3,delta3}],ann2]]

Hence, ann3 is a system of PDEs for the third integral (integration w.r.t. x1, x2, and x3).

Timing[{ann4,delta4}=FindCreativeTelescoping[ann3,Der[x4]];]

(* Test *)
Timing[OreReduce[ann4[[1]]+Der[x4]**delta4[[1,1]],ann3]]

Hence, ann4 is an ODE for the original four-fold integral.

Factor[ann4]

The same differential equation, using the notation (Subscript[θ, x])=x (Subscript[D, x]):

Factor[ChangeOreAlgebra[z**ann4[[1]],OreAlgebra[Euler[z]]]]

Return probability

The following command computes a recurrence for f(n)=p(0)+p(1)+...+p(n), where p(k) is the probability that a random walk ends at the origin after exactly k steps.

Timing[rec4=RecurrenceForReturnProbability[ann4,Der[z],f[n]]]

These are the asymptotic solutions of this recurrence:

Asymptotics[rec4,f[n],Order->3]

The number of excursions in the four-dimensional fcc lattice:

exc4={1,0,24,192,3384,51840};

Timing[data=UnrollRecurrence[rec4,f[n],Table[Sum[exc4[[i+1]]/24^i,{i,0,n}],{n,0,5}],5000];]

Timing[lim4=MyLimit[data,Asymptotics[rec4,f[n],Order->30],n]]

Timing[lim4a=MyLimit[data,Asymptotics[rec4,f[n],Order->32],n]]

The results agree on more than 100 digits, and thus we can conclude that (at least) the first 100 digits are correct:

lim4-lim4a

The return probability therefore is

1-1/lim4a

The five-dimensional face-centred cubic lattice

Differential equation for the LGF

TraditionalForm[HoldForm[Integrate[1/(1-z/10*(Cos[k1]*Cos[k2]+Cos[k1]*Cos[k3]+Cos[k1]*Cos[k4]+Cos[k1]*Cos[k5]+Cos[k2]*Cos[k3]+Cos[k2]*Cos[k4]+Cos[k2]*Cos[k5]+Cos[k3]*Cos[k4]+Cos[k3]*Cos[k5]+Cos[k4]*Cos[k5])),{k1,0,Pi},{k2,0,Pi},{k3,0,Pi},{k4,0,Pi},{k5,0,Pi}]]]

After the substitutions xi = cos (ki) the integrand transforms to :

integrand=1/(1-z/10*(x1*x2+x1*x3+x1*x4+x1*x5+x2*x3+x2*x4+x2*x5+x3*x4+x3*x5+x4*x5))/(Sqrt[1-x1^2]*Sqrt[1-x2^2]*Sqrt[1-x3^2]*Sqrt[1-x4^2]*Sqrt[1-x5^2]);

ann0=Annihilator[integrand,{Der[x1],Der[x2],Der[x3],Der[x4],Der[x5],Der[z]}];

ann0 is a system of PDEs for the integrand.

Timing[{ann1,delta1}=CreativeTelescoping[ann0,Der[x1],{Der[x2],Der[x3],Der[x4],Der[x5],Der[z]}];]

(* Test *)
Timing[OreReduce[MapThread[(#1+Der[x1]**#2)&,{ann1,delta1}],ann0]]

Hence, ann1 is a system of PDEs for the first integral (integration w.r.t. x1).

Timing[{ann2,delta2}=CreativeTelescoping[ann1,Der[x2],{Der[x3],Der[x4],Der[x5],Der[z]}];]

(* Test *)
Timing[OreReduce[MapThread[(#1+Der[x2]**#2)&,{ann2,delta2}],ann1]]

Hence, ann2 is a system of PDEs for the second integral (integration w.r.t. x1 and x2).

To compute the following certificates takes some time, and therefore they are provided in the file fcc5_cert.m.

{ann3,delta3,ann4,delta4,ann5,delta5}=<<"fcc5_cert.m";

Timing[OreReduce[MapThread[(#1+Der[x3]**#2)&,{ann3,delta3}],ann2]]

Hence, ann3 is a system of PDEs for the third integral (integration w.r.t. x1, x2, and x3).

Timing[OreReduce[MapThread[(#1+Der[x4]**#2)&,{ann4,delta4}],ann3]]

Hence, ann4 is a system of PDEs for the fourth integral (integration w.r.t. x1, x2, x3, and x4).

Timing[OreReduce[ann5[[1]]+Der[x5]**delta5[[1]],ann4]]

Hence, ann5 is an ODE for the original five-fold integral.

Factor[ann5]

The same differential equation, using the notation (Subscript[θ, x])=x (Subscript[D, x]):

Factor[ChangeOreAlgebra[z^2**ann5[[1]],OreAlgebra[Euler[z]]]]

Return probability

The following command computes a recurrence for f(n)=p(0)+p(1)+...+p(n), where p(k) is the probability that a random walk ends at the origin after exactly k steps.

Timing[rec5=RecurrenceForReturnProbability[ann5,Der[z],f[n]]]

These are the asymptotic solutions of this recurrence:

Asymptotics[rec5,f[n],Order->3]

The number of excursions in the five-dimensional fcc lattice:

exc5={1,0,40,480,11880,281280,7506400,210268800};

Timing[data=UnrollRecurrence[rec5,f[n],Table[Sum[exc5[[i+1]]/40^i,{i,0,n}],{n,0,7}],5000];]

Timing[lim5=MyLimit[data,Asymptotics[rec5,f[n],Order->34],n]]

Timing[lim5a=MyLimit[data,Asymptotics[rec5,f[n],Order->36],n]]

The results agree on more than 100 digits, and thus we can conclude that (at least) the first 100 digits are correct:

lim5-lim5a

The return probability therefore is

1-1/lim5a

The six-dimensional face-centred cubic lattice

Differential equation for the LGF

combs=Total[Times@@@Subsets[Symbol["x"<>ToString[#]]&/@Range[6],{2}]]

integrand=1/(1-z/Length[combs]*combs)/Product[Sqrt[1-Symbol["x"<>ToString[i]]^2],{i,6}]

ann0=Annihilator[integrand,{Der[x1],Der[x2],Der[x3],Der[x4],Der[x5],Der[x6],Der[z]}];

ann0 is a system of PDEs for the integrand.

Timing[{ann1,delta1}=CreativeTelescoping[ann0,Der[x1],{Der[x2],Der[x3],Der[x4],Der[x5],Der[x6],Der[z]}];]

(* Test *)
Timing[OreReduce[MapThread[(#1+Der[x1]**#2)&,{ann1,delta1}],ann0]]

Hence, ann1 is a system of PDEs for the first integral (integration w.r.t. x1).

Timing[{ann2,delta2}=CreativeTelescoping[ann1,Der[x2],{Der[x3],Der[x4],Der[x5],Der[x6],Der[z]}];]

(* Test *)
Timing[OreReduce[MapThread[(#1+Der[x2]**#2)&,{ann2,delta2}],ann1]]

Hence, ann2 is a system of PDEs for the second integral (integration w.r.t. x1 and x2).

{ann3,delta3,ann4,delta4,ann5,delta5,ann6,delta6}=<<"fcc6_cert.m";

Timing[OreReduce[MapThread[(#1+Der[x3]**#2)&,{ann3,delta3}],ann2]]

(* A quicker (but non-rigorous) correctness check. *)
subs={x4->92,x5->55,x6->-74,z->23};
Timing[{ann3a,delta3a}=OrePolynomialSubstitute[#,subs]&/@{ann3,delta3};
OreReduce[MapThread[(#1+Der[x3]**#2)&,{ann3a,delta3a}],ann2,OrePolynomialSubstitute->subs]]

Hence, ann3 is a system of PDEs for the third integral (integration w.r.t. x1, x2, and x3).

Timing[OreReduce[MapThread[(#1+Der[x4]**#2)&,{ann4,delta4}],ann3]]

(* A quicker (but non-rigorous) correctness check. *)
subs={x5->55,x6->-74,z->23};
Timing[{ann4a,delta4a}=OrePolynomialSubstitute[#,subs]&/@{ann4,delta4};
OreReduce[MapThread[(#1+Der[x4]**#2)&,{ann4a,delta4a}],ann3,OrePolynomialSubstitute->subs]]

Hence, ann4 is a system of PDEs for the fourth integral (integration w.r.t. x1, x2, x3, and x4).

Timing[OreReduce[MapThread[(#1+Der[x5]**#2)&,{ann5,delta5}],ann4]]

(* A quicker (but non-rigorous) correctness check. *)
subs={x6->-74,z->23};
Timing[{ann5a,delta5a}=OrePolynomialSubstitute[#,subs]&/@{ann5,delta5};
OreReduce[MapThread[(#1+Der[x5]**#2)&,{ann5a,delta5a}],ann4,OrePolynomialSubstitute->subs]]

Hence, ann5 is a system of PDEs for the fifth integral (integration w.r.t. x1, x2, x3, x4, and x5).

Timing[OreReduce[ann6[[1]]+Der[x6]**delta6[[1]],ann5]]

(* A quicker (but non-rigorous) correctness check. *)
subs={z->23};
Timing[{ann6a,delta6a}=OrePolynomialSubstitute[#,subs]&/@{ann6,delta6};
OreReduce[MapThread[(#1+Der[x6]**#2)&,{ann6a,delta6a}],ann5,OrePolynomialSubstitute->subs]]

Hence, ann6 is an ODE for the original six-fold integral.

Factor[ann6]

Return probability

The following command computes a recurrence for f(n)=p(0)+p(1)+...+p(n), where p(k) is the probability that a random walk ends at the origin after exactly k steps.

Timing[rec6=RecurrenceForReturnProbability[ann6,Der[z],f[n]];]

These are the asymptotic solutions of this recurrence:

Asymptotics[rec6,f[n],Order->3]

The number of excursions in the six-dimensional fcc lattice:

exc6={1,0,60,960,30780,996480,36560400,1430553600,59089923900,2543035488000,113129280527760,5170796720812800,241741903350301200};

Timing[data=UnrollRecurrence[rec6,f[n],Table[Sum[exc6[[i+1]]/60^i,{i,0,n}],{n,0,12}],5000];]

Timing[lim6=MyLimit[data,Asymptotics[rec6,f[n],Order->34],n]]

Timing[lim6a=MyLimit[data,Asymptotics[rec6,f[n],Order->36],n]]

The results agree on more than 100 digits, and thus we can conclude that (at least) the first 100 digits are correct:

lim6-lim6a

The return probability therefore is

1-1/lim6a