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

NN=4;
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,z,w,α,β];

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

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

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

Timing[ann1=FindCreativeTelescoping[ann0,Der[x1]][[1]];]

Timing[ann2=FindCreativeTelescoping[ann1,Der[x2]][[1]];]

Timing[ann3=FindCreativeTelescoping[ann2,Der[x3]][[1]];]

Timing[ann4=FindCreativeTelescoping[ann3,Der[x4]][[1]];]

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

ann4=ToExpression[Import[NotebookDirectory[]<>"Data-N4M3-Integral.txt"]];

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

ODEDiv2=ann4[[1]];

Compute the ODE for R(z).
ODEinD - in terms of the derivation operator D
ODEinTheta - in terms of the derivation operator θ - Order 8, Degree 32 (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 4.2!  (To be displayed at the end of this notebook)
Order 8, Degree 16

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 16

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=20;
ClearAll[a];

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

SeriesIniSol=Solve[Join[Table[Coefficient[SeriesIni,w,i]==0,{i,1,MAX}],{a[0]==1}],Table[a[i],{i,0,MAX}]]

SeqListIni=Table[SeriesIniSol[[1,k,2]],{k,1,Length@SeriesIniSol[[1]]}]

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

The above values of (OverTilde[r])(n) can be compared with those computed directly.

SeqListIniComputation={};

MAX=20;

For[n=0,n≤MAX,n++,
coord=Select[Tuples[Table[i,{i,0,n}],NN],Total[#]==n&];
size=Length@coord;
p=Sum[Multinomial[Sequence@@(2coord[[i]])]*Product[Binomial[2n-2coord[[i,j]],n-coord[[i,j]]],{j,1,NN}],{i,1,size}];
AppendTo[SeqListIniComputation,p];
];

SeqListIniComputation

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

Bound=5000;

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

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

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

ClearAll[Diff];
ODEGuess=GuessMinDE[Take[SeqList,300],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]]

Compute the asymptotics for (OverTilde[r])(n).

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

Ind=Reverse[Table[Floor[Bound/i],{i,1,3}]]
Table[N[((seq[Ind[[i]]])/(AsyList[[1]]/.{α->Ind[[i]]}))],{i,1,Length@Ind}]
Table[N[((seq[Ind[[i]]])/(AsyList[[2]]/.{α->Ind[[i]]}))],{i,1,Length@Ind}]
Table[N[((seq[Ind[[i]]])/(AsyList[[3]]/.{α->Ind[[i]]}))],{i,1,Length@Ind}]
Table[N[((seq[Ind[[i]]])/(AsyList[[4]]/.{α->Ind[[i]]}))],{i,1,Length@Ind}]
Table[N[((seq[Ind[[i]]])/(AsyList[[11]]/.{α->Ind[[i]]}))],{i,1,Length@Ind}]
Table[N[((seq[Ind[[i]]])/(AsyList[[12]]/.{α->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 for (OverTilde[R])(z) in Theorem 4.2

ODENormalizedinTheta/.{w->z}