Multi-headed Lattice Green Function (N = 5, M = 3)

NN=5;
MM=3;

Recall some basic definitions in the paper:

(Subscript[P, M,N])(z):=(1/(((2 π))^N))(∫_(-π)^π)...(∫_(-π)^π)(1/(1 - (z/((({{N}, {M}})))) (Subscript[σ, M])(cos (Subscript[θ, 1]), ..., cos (Subscript[θ, N]))))d(Subscript[θ, 1])...d(Subscript[θ, N])

(Subscript[R, M,N])(z):=(Subscript[P, M,N])((2^M)(({{N}, {M}}))z)   and  (Subscript[R, M,N])(z)=(∑_(n≥0))(Subscript[r, M,N])(n)(z^n)

Also, for M odd or M=N, we always have r(2n+1)=0. Hence, define
(Subscript[(OverTilde[r]), M,N])(n):=(Subscript[r, M,N])(2n)  and  (Subscript[(OverTilde[R]), M,N])(z):=(∑_(n≥0))(Subscript[(OverTilde[r]), M,N])(n)(z^n)=(∑_(n≥0))(Subscript[r, M,N])(2n)(z^n) 

Our goal is to find:

Case 1. M even and M≠N:
- recurrences (REC) for r(n) or differential equations (ODE) for R(z).

Case 2. M odd or M=N:
- recurrences (REC) for (OverTilde[r])(n) or differential equations (ODE) for (OverTilde[R])(z).

Command: UnrollRecurrence
Generate a sequence from recurrence & initial values (Koutschan’s implementation).

(* 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];
];

Load RISC packages.

<<RISC`HolonomicFunctions`
<<RISC`Asymptotics`
<<RISC`Guess`

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Asymptotics Package version 0.3
written by Manuel Kauers
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Package GeneratingFunctions version 0.9 written by Christian Mallinger
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Guess Package version 0.52
written by Manuel Kauers
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Apply creative telescoping to R(z/(2^M)).

ClearAll[x1,x2,x3,x4,x5,z,w,α,β];

SymmetricPolynomial[3,{x1,x2,x3,x4,x5}]

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

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

Timing[{ann1,delta1}=CreativeTelescoping[ann0,Der[x1]];]

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

Alternatively, you may import the value of ann2 from an external file.

{ann2,delta2}=ToOrePolynomial[ToExpression[Import[NotebookDirectory[]<>"Data-N5M3-Integral-ann2.txt"]]];

Now you need to import the annihilating operators from external files. In particular, the certification of the telescopers has been done on a server; here we only display the certification with specific substitutions of the variables.

{ann3,delta3}=ToOrePolynomial[ToExpression[Import[NotebookDirectory[]<>"Data-N5M3-Integral-ann3.txt"]]];

(*Certify the telescopers for the third integral numerically (and hence in a nonrigorous way): Output needs to be {0, 0, 0, 0}*)
subs={x4->-66,x5->497,z->333};
{ann3a,delta3a}=OrePolynomialSubstitute[#, subs]&/@{ann3,delta3};
Timing[OreReduce[MapThread[(#1+Der[x3]**#2)&,{ann3a,delta3a}],ann2,OrePolynomialSubstitute->subs]]

{ann4,delta4}=ToOrePolynomial[ToExpression[Import[NotebookDirectory[]<>"Data-N5M3-Integral-ann4.txt"]]];

(*Certify the telescopers for the fourth integral numerically (and hence in a nonrigorous way): Output needs to be {0, 0, 0, 0}*)
subs={x5->-129,z->654};
{ann4a,delta4a}=OrePolynomialSubstitute[#, subs]&/@{ann4,delta4};
Timing[OreReduce[MapThread[(#1+Der[x4]**#2)&,{ann4a,delta4a}],ann3,OrePolynomialSubstitute->subs]]

{ann5,delta5}=ToOrePolynomial[ToExpression[Import[NotebookDirectory[]<>"Data-N5M3-Integral-ann5.txt"]]];

(*Certify the telescopers for the fifth integral numerically (and hence in a nonrigorous way): Output needs to be {0}*)
subs={z->11};
{ann5a,delta5a}=OrePolynomialSubstitute[#, subs]&/@{ann5,delta5};
Timing[OreReduce[MapThread[(#1+Der[x5]**#2)&,{ann5a,delta5a}],ann4,OrePolynomialSubstitute->subs]]

ann5 gives an ODE for R(z/(2^M)).

ODEDiv2=ann5[[1]];

Compute the ODE for R(z).
ODEinD - in terms of the derivation operator D
ODEinTheta - in terms of the derivation operator θ - Order 14, Degree 110 (Refer to Table 1)

ODETemp=NormalizeCoefficients[DFiniteSubstitute[{ODEDiv2},{z->w*(2^MM)},Algebra->OreAlgebra[Der[w]]][[1]]];

ODEinD=NormalizeCoefficients[DFiniteSubstitute[{ODETemp},{w->z},Algebra->OreAlgebra[Der[z]]][[1]]];

ODEinTheta=NormalizeCoefficients[ChangeOreAlgebra[z**ODEinD,OreAlgebra[Euler[z]]]];

ODEinThetaOrder=OrePolynomialDegree[ODEinTheta,Euler[z]]

ODEinThetaDegree=Max[Exponent[OrePolynomialListCoefficients[ODEinTheta],z]]

Since M=3 is odd, we move on to the ODE for (OverTilde[R])(z)=R((z^(1/2))).
ODENormalizedinTheta gives the ODE in Theorem 5.1!  (To be displayed at the end of this notebook)
Order 14, Degree 55

ODENormalizedinD=NormalizeCoefficients[DFiniteSubstitute[{ODEinD},{z->(w^(1/2))},Algebra->OreAlgebra[Der[w]]][[1]]];

ODENormalizedinTheta=NormalizeCoefficients[ChangeOreAlgebra[w**ODENormalizedinD,OreAlgebra[Euler[w]]]];

ODENormalizedinThetaOrder=OrePolynomialDegree[ODENormalizedinTheta,Euler[w]]

ODENormalizedinThetaDegree=Max[Exponent[OrePolynomialListCoefficients[ODENormalizedinTheta],w]]

Get the REC for (OverTilde[r])(n).
Order 55

RECNormalizedinS=NormalizeCoefficients[DFiniteDE2RE[{ODENormalizedinD},{w},{α}][[1]]];

RecNormalizedinSOrder=OrePolynomialDegree[RECNormalizedinS,S[α]]

We may also write this REC explicitly.

ClearAll[Seq];
SeqNormalized=ApplyOreOperator[RECNormalizedinS,Seq[α]];

The initial values of (OverTilde[r])(n) are also produced by the ODE for (OverTilde[R])(z).

MAX=ODENormalizedinThetaDegree;
ClearAll[a];

SeriesIni=ApplyOreOperator[ODENormalizedinTheta,Sum[a[n](w^n),{n,0,MAX}]];

SeqListIni={1,80};

For[k=2,k≤MAX,k++,
{
eqn=Coefficient[SeriesIni,w,k]==0/.Table[a[i]->SeqListIni[[i+1]],{i,0,k-1}];
eqnsol=Solve[eqn,a[k]][[1,1,2]];
AppendTo[SeqListIni,eqnsol];
}
];

SeqListIni

seq[n_]:=SeqListIni[[n+1]];

seq[54]

Alternatively, you may import these initial values from an external file.

SeqListIni=ToExpression[Import[NotebookDirectory[]<>"Data-N5M3-Integral-Initial-Values.txt"]];

seq[n_]:=SeqListIni[[n+1]];

Let us the generate a list of (OverTilde[r])(n).

Bound=10000;

SeqList=UnrollRecurrence[SeqNormalized,Seq[α],SeqListIni,Bound];

seq[n_]:=SeqList[[n+1]];

Guess a Minimal ODE for (OverTilde[R])(z).
Its order is 14, and is identical to that of the ODE in Theorem 5.1 (ODENormalizedinTheta).

ClearAll[Diff];
ODEGuess=GuessMinDE[Take[SeqList,800],Diff[z]];
ODEGuessinD=NormalizeCoefficients[ToOrePolynomial[ODEGuess/.{Derivative[k_][Diff][z]->((Der[z])^k)}/.{Diff[z]->1}]];

ODEGuessinTheta=NormalizeCoefficients[ChangeOreAlgebra[z**ODEGuessinD,OreAlgebra[Euler[z]]]];

ODEGuessinThetaOrder=OrePolynomialDegree[ODEGuessinTheta,Euler[z]]

Guess a Minimal REC for (OverTilde[r])(n).
SeqfromRECGuess gives the REC in Theorem 5.2!  (To be displayed at the end of this notebook)
REC: Order 8
ODE: Order 69, Degree 8

RECGuess=GuessMinRE[Take[SeqList,800],Seq[α]];
RECGuessinS=NormalizeCoefficients[ToOrePolynomial[RECGuess/.{Seq[k_]->((S[α])^(k-α))}]];

RECGuessinSOrder=OrePolynomialDegree[RECGuessinS,S[α]]

ODEfromRECGuessinD=NormalizeCoefficients[DFiniteRE2DE[{RECGuessinS},{α},{z}][[1]]];

ODEfromRECGuessinTheta=NormalizeCoefficients[ChangeOreAlgebra[z**ODEfromRECGuessinD,OreAlgebra[Euler[z]]]];

ODEfromRECGuessinThetaOrder=OrePolynomialDegree[ODEfromRECGuessinTheta,Euler[z]]

ODEfromRECGuessinThetaDegree=Max[Exponent[OrePolynomialListCoefficients[ODEfromRECGuessinTheta],z]]

We may also write this REC explicitly.

ClearAll[Seq];
SeqfromRECGuess=ApplyOreOperator[RECGuessinS,Seq[α]];

SeqfromRECGuessList=UnrollRecurrence[SeqfromRECGuess,Seq[α],Take[SeqList,RECGuessinSOrder],200];

Prove the minimal REC for (OverTilde[r])(n).

RECCompare=DFinitePlus[{RECNormalizedinS},{RECGuessinS}][[1]];

Compute the largest positive integral root of the leading coefficient in the recurrence RECCompare.

LeadCoeff=RECCompare[[1,1,1]];
LeadCoeffRoot=Solve[LeadCoeff==0,α][[All,1,2]]

There are no positive integral roots in our case.

Select[Select[LeadCoeffRoot,IntegerQ],#>0&]

RECCompareOrder=LeadingExponent[RECCompare][[1]]

CheckNum=RECCompareOrder+20;
Take[SeqList,CheckNum]-Take[SeqfromRECGuessList,CheckNum]

Compute the asymptotics for (OverTilde[r])(n).
This code requires too much memory so we have run it on a server. The output indicates the following asymptotic formula.

AsyList=Asymptotics[SeqfromRECGuess,Seq[α]];
N[AsyList]

Ind=Reverse[Table[Floor[Bound/i],{i,1,3}]]
Table[N[((seq[Ind[[i]]])/(((6400^α)(α^(-5/2)))/.{α->Ind[[i]]}))],{i,1,Length@Ind}]

Approximate the Polya number.

AtOne=N[Sum[seq[n]*((((1/((2^MM)Binomial[NN,MM]))))^(2n)),{n,0,Bound}],11]
N[1-(1/AtOne),10]

Display the ODE in Theorem 5.1

-ODENormalizedinTheta/.{w->z}

Display the REC in Theorem 5.2

Collect[Expand[SeqfromRECGuess],Seq[_]]