**Staircase polygons, a simplified model for self-avoiding walks**

*Frederic Chyzak*

(Version of January 16, 1997)

A combinatorial problem in statistical physics is that of counting
*self-avoiding walks*
in an integer lattice. In particular, the cases of dimensions 2 and 3 are related to problems of phase transitions in the classical Ising model of statistical physics. Solving the general model for self-avoiding walks is an old open problem so that one focusses on approximate models. Exactly solvable models exist for various models obtained by constraining the admissible walks. More and more is known about self-avoiding walks by relaxing these constraints.

Here we follow a thread started by A. J. Guttmann and T. Prellberg [
*Staircase polygons, elliptic integrals, Heun functions, and lattice Green functions*
, (1993), Physical Review E, (47) 4, R2233-2236]. We show how the
__Gfun__
package may be used to derive in a matter of seconds differential equations that were previously guessed in the article above, and how to establish rigorously the singular structure of the ODE's. As a by-product, we derive effective asymptotic estimates. Early results on this subject date back to G. Polya [
*On the number of certain lattice polygons*
, (1969), Journal of Combinatorial Theory, Series A, (6), 102-105].

**Introduction**

A strong but simple constraint on the walks is to require them to proceed by positive steps. Throughout the remainder of this session, we only consider such paths with positive steps.

A
*-dimensional staircase polygon*
is a special self-avoiding walk in the
-dimensional integer lattice where all positive steps precede all negative steps. This can be interpreted as an (unordered) pair of non-crossing paths of length
from the same origin to the same end. Here is an example of length
in dimension
:

`> `
**plot([[0,0],[2,0],[2,1],[4,1],[4,3],[2,3],[2,2],[1,2],[1,1],[0,1],[0,0]]);**

A
*-dimensional festoon*
is obtained from a sequence of staircase polygons, with the possible introduction of unit steps.

Here is an example in dimension 2, involving 3 polygons (thick lines) and 5 unit steps (thin lines).

`> `
**plots[display](
plot([[0,0],[2,0],[2,2],[1,2],[1,1],[0,1],[0,0]],thickness=3),
plots[textplot]([1/2,3/2,`+`],font=[TIMES,BOLD,18]),
plots[textplot]([5/2,1/2,`-`],font=[TIMES,BOLD,18]),
plot([[2,2],[3,2]]),
plot([[3,2],[5,2],[5,5],[6,5],[6,6],[4,6],[4,4],[3,4],[3,2]],thickness=3),
plots[textplot]([7/2,9/2,`+`],font=[TIMES,BOLD,18]),
plots[textplot]([11/2,7/2,`-`],font=[TIMES,BOLD,18]),
plot([[6,6],[7,6],[7,7]]),
plot([[7,7],[10,7],[10,8],[11,8],[11,9],[8,9],[8,8],[7,8],[7,7]],thickness=3),
plots[textplot]([21/2,15/2,`+`],font=[TIMES,BOLD,18]),
plots[textplot]([15/2,17/2,`-`],font=[TIMES,BOLD,18]),
plot([[11,9],[13,9]]),
scaling=constrained,title=`A 2-dimensional festoon`);**

We orient each of the staircase polygon (i.e., we distinguish between its two consisting paths), so that a festoon becomes a sequence of oriented polygons, unit steps remaining non-oriented. On the previous display, this orientation is marked by and signs. In this way, we interpret a festoon as an (oriented) pair of paths of length with same origin and end.

Denote by the number of staircase polygons in the -dimensional integer lattice, and call the ordinary generating function

.

This function has a simple relation to the ordinary generating function

of the number of (ordered) pairs of paths of length with same origin and end. More precisely,

,

since any pair of paths can be viewed as a sequence of pairs of non-crossing paths. Equivalently,

.

In this session, we proceed to compute the singularities and asymptotics of for small dimensions.

**Dimension 2: an explicit case**

The case of dimension 2 is simple and yields explicit closed forms for generating functions and counting sequences, already known to Polya. It is usually solved by hand using a Vandermonde convolution. It is readily solved by Maple. The number of pairs of path of length with same origin and end is:

`> `
**p[2,n]=Sum(binomial(n,k)^2,k=0..n);**

because a path is determined by the choice of its steps in one of the directions of the lattice. The number of pairs of paths is thus given by the central binomial coefficient:

`> `
**p[2,n]=normal(convert(sum(binomial(n,k)^2,k=0..n),factorial),expanded);**

The generating function of these numbers is well-known:

`> `
**P[2](z)=sum(op(2,")*z^n,n=0..infinity);**

Another closed form for the coefficients of this series are gotten by Newton's binomial expansion:

`> `
**p[2,n]=binomial(-1/2,n)*(-4)^n;**

The following asymptotics for is obtained by Maple:

`> `
**asympt_p[2]:=expand(simplify(subs(cos(Pi*n)=(-1)^n,convert(asympt(op(2,"),n,2),polynom)),symbolic));**

From the value found for , we get the generating function for the non-crossing paths:

`> `
**Q[2](z)=(1+2*z-1/op(2,"""))/2;**

which Maple easily expands into:

`> `
**series(op(2,"),z,10);**

A closed form for the coefficients of this series are gotten by Newton's binomial expansion:

`> `
**q[2,n]=delta[n,0]/2+delta[n,1]-binomial(1/2,n)*(-4)^n/2;**

where denotes the Kronecker symbol, if and only if , 0 otherwise. The following evaluation corroborates this expansion:

`> `
**seq(-binomial(1/2,i)*(-4)^i/2,i=2..9);**

The following asymptotics for is obtained by Maple:

`> `
**asympt_q[2]:=expand(simplify(subs(cos(Pi*n)=(-1)^n,convert(asympt(-(-1)^n*4^n*binomial(1/2,n)/2,n,2),polynom)),symbolic));**

We check this asymptotic estimate using
__Gfun__
.

`> `
**with(gfun);**

The following first order differential equation is satisfied by the generating function :

`> `
**2*f(z)+(1-4*z)*diff(f(z),z);**

`> `
**normal(eval(subs(f(z)=sqrt(1-4*z),")));**

We compute a recursion on the coefficients
of the series using
__gfun[diffeqtorec]__
.

`> `
**diffeqtorec("",f(z),u(n));**

Next,
__gfun[rectoproc]__
returns a procedure that computes the coefficients in an efficient way.

`> `
**coeff_of_Q:=rectoproc({",u(0)=1,u(1)=1,u(2)=1},u(n));**

We finally check our asymptotic estimate: the following ratio goes to 1 when goes to infinity.

`> `
**for i from 100 to 1000 by 100 do i=evalf(coeff_of_Q(i)/subs(n=i,asympt_q[2]),10) od;**

Using
__series__
to expand
and compute the coefficients
would be less efficient and require more time and space.

In higher dimensions, the problem no longer has a known explicit solution. A natural tool to attack it are
*Bessel generating functions*
. The Bessel generating function of a sequence
of numbers is defined as the sum

.

For instance, the Bessel generating function of the constant sequence
is given by the
__modified Bessel function__

,

which we denote by j(z). The Bessel generating function of the numbers of pairs of paths with same origin and end and positive steps in a -dimensional lattice is given by the product

,

where marks the steps in dimension . It follows that the Bessel generating function of the numbers with respect to the length of the paths is .

**Higher dimensional cases: computing differential equations**

In dimension 2, the generating function
is algebraic. In higher dimensions,
belongs to the larger class of
*holonomic functions*
. A function
is called holonomic when it satisfies a linear differential equation with rational function coefficients. This class of functions benefits from numerous closure properties for which algorithms have been implemented in the package
__Gfun__
. In particular, this class is closed under sum, product, Borel and inverse Borel transforms. The
__Borel transform__
defined in Gfun is related to the
__inverse Laplace transform__
: formally applied on a formal power series, it is

.

This is enough for us to compute differential equations satisfied by the .

**Dimension 2 revisited**

This case was solved above in an explicit way.

We proceed to get , and more generally , by the formula

.

We first compute a differential equation satisfied by starting from a recurrence satisfied by its coefficients :

`> `
**diff_eq_j:=rectodiffeq({(n+1)^2*u(n+1)=u(n),u(0)=1},u(n),j(z));**

We compute equations satisfied by the powers
using
__gfun[poltodiffeq]__
.

`> `
**power_j:=proc(d) option remember; poltodiffeq(j(z)^d,[diff_eq_j],[j(z)],j(z)) end:**

For instance, in the case , the following system is satisfied by :

`> `
**power_j(2);**

We now derive a differential system satisfied by
, the ordinary generating function of the
, by computing two inverse Borel transforms with
__gfun[invborel]__
.

`> `
**double_inv_borel:=proc(sys,f,z) option remember; invborel(invborel(sys,f(z),diffeq),f(z),diffeq) end:**

In the case , we have:

`> `
**double_inv_borel(power_j(2),j,z);**

The equation is solve by
__dsolve__
.

`> `
**dsolve(",j(z));**

The possible singularities of a holonomic equation are given by its leading coefficient (on the previous equation,
) and these singularities yield the asymptotic form of coefficients by the
*method of singularity analysis*
. The singularity read on the previous equation is in
, so that the exponential behaviour of
is
, which we found in the previous section.

In the case of dimension 2, we read the singularity by:

`> `
**solve(coeff(op(remove(type,"",equation)),diff(j(z),z)),z);**

**Dimension 3**

We perform the same calculations for dimensions 3 and higher.

System satisfied by .

`> `
**power_j(3);**

System satisfied by .

`> `
**double_inv_borel(",j,z);**

Order and possible singularities.

`> `
**op(sort([solve(coeff(op(remove(type,",equation)),diff(j(z),z,z)),z)]));**

**Dimension 4**

System satisfied by .

`> `
**power_j(4);**

System satisfied by .

`> `
**double_inv_borel(",j,z);**

Order and possible singularities.

`> `
**op(sort([solve(coeff(op(remove(type,",equation)),diff(j(z),z$3)),z)]));**

**Dimension 5**

System satisfied by .

`> `
**power_j(5);**

System satisfied by .

`> `
**double_inv_borel(",j,z);**

Order and possible singularities.

`> `
**op(sort([solve(coeff(op(remove(type,",equation)),diff(j(z),z$4)),z)]));**

**Dimension 6**

System satisfied by .

`> `
**power_j(6);**

System satisfied by .

`> `
**double_inv_borel(",j,z);**

Order and possible singularities.

`> `
**op(sort([solve(coeff(op(remove(type,",equation)),diff(j(z),z$5)),z)]));**

**Dimension 7**

System satisfied by

`> `
**power_j(7);**

System satisfied by .

`> `
**double_inv_borel(",j,z);**

Order and possible singularities.

`> `
**op(sort([solve(coeff(op(remove(type,",equation)),diff(j(z),z$6)),z)]));**

These experimental results suggest a pattern for the singularities of the equation in dimension . More specifically, the leading coefficient of the differential equation is given by

when , and by

when .

**The singular structure of the ODE**

The equation used to describe introduces a "parasitic" function, , also expressed in terms of a modified Bessel function.

`> `
**op(remove(type,diff_eq_j,equation));**

Both and are solutions found my Maple:

`> `
**dsolve(",j(z));**

The equation above is satisfied by
*any*
linear combination of
and
. In the same way, the equation used to described
is satisfied by
*any*
linear combination of the products
. For instance, this is easily checked for the cases
and
.

**Dimension 2**

Here is the differential equation satisfied by :

`> `
**op(remove(type,power_j(2),equation));**

Maple is no longer able to solve in terms of special functions:

`> `
**dsolve(",j(z));**

(The development version of Maple is able to solve.)

We successively plug , and into the previous equation:

`> `
**op(remove(type,power_j(2),equation)):**

`> `
**normal(eval(subs(j=unapply(BesselI(0,2*sqrt(z))^2,z),")));**

`> `
**normal(eval(subs(j=unapply(BesselK(0,2*sqrt(z))^2,z),"")));**

`> `
**normal(eval(subs(j=unapply(BesselI(0,2*sqrt(z))*BesselK(0,2*sqrt(z)),z),""")));**

Each evaluation yields 0, proving that the function is a solution.

**Dimension 3**

Here is the differential equation satisfied by :

`> `
**op(remove(type,power_j(3),equation));**

We successively plug , , and into the previous equation:

`> `
**normal(eval(subs(j=unapply(BesselI(0,2*sqrt(z))^3,z),")));**

`> `
**normal(eval(subs(j=unapply(BesselI(0,2*sqrt(z))^2*BesselK(0,2*sqrt(z)),z),"")));**

`> `
**""":**

`> `
**normal(eval(subs(j=unapply(BesselI(0,2*sqrt(z))*BesselK(0,2*sqrt(z))^2,z),")));**

`> `
**normal(eval(subs(j=unapply(BesselK(0,2*sqrt(z))^3,z),"")));**

The ordinary generating function is the double inverse Borel transform of . In the same way as above, the differential equation that vanishes at also vanishes at each of the double inverse Borel transforms . Those transform have a nice integral representation, namely

.

This integral is always convergent at its lower bound. At infinity however, the integrand behaves like

.

It follows that either and the integral is defined for any non-negative , or and the integral is defined for between 0 and . We thus obtain all positive singularities of the conjecture.

**Theorem:**
Any homogeneous linear ODE satisfied by the generating function
* of festoons on the *
*-dimensional integer lattice is singular at each *
*, for *
. More specifically, the function
* is singular at *
*.*

**Computing the asymptotics: a connection problem**

The function
is analytic in a neighbourhoud of 0, where it is given by its Taylor expansion. On the other hand, its dominant singularity (i.e., its singularity closest to 0) and its asymptotic behaviour there are a priori not known. It is natural to estimate the asymptotics of
at the smallest positive singularity,
, of the equation defining this function. To do so, we use
*singularity analysis*
.
First, we evaluate the Taylor series at 0.
Next, we algorithmically compute a basis of all possible asymptotic behaviours at
of solutions of the equation.
Finally, we numerically connect the Taylor series at 0 and a generic linear combination of these asymptotics at
, i.e., we equate values at a point between 0 and
.

A common case of singularity of holonomic function is that of a
*regular singular point*
. This is a singular point where the function has a polynomial behaviour: at such points, the theory expects an expansion of the form
, where
is a formal power series. It is therefore natural to look for such an asymptotics by a method of "undetermined exponent". The possible values for
are given by a polynomial, the
*indicial polynomial*
: a zero
of multiplicity
of the indicial polynomial indicates solutions with the behaviours
,
, ..., and
. The following procedure computes the indicial polynomial of a differential equation.

`> `
**indicial_poly:=proc(sys,f,z) local zero; global eta; option remember;
zero:=collect(normal(eval(subs(f=proc(z) global eta; z^eta end,op(remove(type,sys,equation)))))/z^eta,z,normal);
factor(primpart(coeff(zero,z,ldegree(zero,z)),eta))
end:**

As a first remark, note that is the only solution of the equation analytic at 0 and such that . To obtain this, let us compute the indicial polynomials at 0 for the first dimensions:

`> `
**for i from 2 to 7 do i=indicial_poly(double_inv_borel(power_j(i),j,z),j,z) od;**

Therefore, the equation satisfied by in dimension has one solution which is analytic at 0 and that are singular at 0, with logarithmic singularities.

**Connection and asymptotics in dimension 3**

The smallest singularity is .

`> `
**op(sort([solve(coeff(op(remove(type,double_inv_borel(power_j(3),j,z),equation)),diff(j(z),z,z)),z)]));**

We first expand the series at 0 to evaluate at, say, . Next, we recenter the equation around to study the local behaviour of the solutions there more explicitly.

To evaluate the Taylor series, we compute a truncation of it. We actually need truncations to several orders to illustrate the slow convergence of the series close to its singularity. We base on the following identity to compute such truncations:

,

where denotes the coefficient of in the series and is any series .

We start from the equation for :

`> `
**double_inv_borel(power_j(3),j,z);**

We compute :

`> `
**algebraicsubs(",j=z*t,j(z));**

We divide by :

`> `
**poltodiffeq(j(z)/(1-z),["],[j(z)],j(z));**

We extract its coefficients:

`> `
**rec_j:=diffeqtorec(",j(z),u(n));**

And we compute a procedure, which to associates the truncated series at to the order .

`> `
**series_j:=rectoproc(rec_j,u(n),params=[t]):**

Next, the series can be computed formally at :

`> `
**normal(series_j(4,t));**

or numerically at a rational number:

`> `
**series_j(4,1/16);**

(This is more efficient than computing the series, and next evaluating it numerically.)

The following plots display the Taylor expansion for increasing an number of terms. Note the behaviour of the series close to the singularity when we increase the number of terms.

`> `
**pl[Taylor]:=plots[display](
plot('series_j'(10,z),z=0..1/9,y=0..5,style=point,color=red),
plot('series_j'(20,z),z=0..1/9,y=0..5,style=point,color=green),
plot('series_j'(100,z),z=0..1/9,y=0..5,style=point,color=blue),
title=`Taylor series (red=10 terms, green=20 terms, blue=100 terms)`):**

`> `
**pl[Taylor];**

**Computing a basis of asymptotics expansions**

To compute asymptotic expansions at the singularity , we find it convenient to center the equation around . We thus make the change of variable , i.e., (and ). To perform this rational change of variable, we appeal to the closure of the class of holonomic functions under algebraic substitution:

`> `
**d_eq:=collect(subs({j=f,z=u},algebraicsubs(double_inv_borel(power_j(3),j,z),j-(1-z)/9,j(z))),diff,factor);**

Let us compute the indicial polynomial at the singularity (now in 0):

`> `
**indicial_poly({d_eq},f,u);**

Hence, we expect one regular solution at 0, and one singular with logarithmic behaviour.

We therefore expect an analytic solution and a logarithmic solution. This is confirmed by the following formal expansion:

`> `
**Order:=3: dsolve(d_eq,f(u),series);**

We compute but do not display the series for a larger order, so as to increase the precision in the numerical calculations below.

`> `
**Order:=15: dsol_ser_f:=dsolve(d_eq,f(u),series):**

Note that there exist algorithms to perform the previous calculation much more efficiently.

The following plot displays the asymptotic behaviours at of the solutions.

`> `
**plots[display](
plot(subs(u=1-9*z,convert(subs({_C1=1,_C2=0},op(2,dsol_ser_f)),polynom)),
z=0..1/5,y=0..5,color=blue),
plot(subs(u=1-9*z,convert(subs({_C1=0,_C2=-1},op(2,dsol_ser_f)),polynom)),
z=0..1/9,y=-2..5,color=green),
plot([[1/9,-2],[1/9,5]],color=black),
title=`Asymptotic behaviours at 1/9 (blue=analytic, green=logarithmic)`);**

The function is a linear combination of both behaviours.

To perform the connection, we determine values for and above so that and agree when takes values around .

`> `
**Digits:=25:**

We graphically determine .

`> `
**plot(subs(_C1=x,map(solve,{seq(series_j(50,(1-(.5+.01*i))/9)=eval(subs(u=.5+.01*i,convert(op(2,dsol_ser_f),polynom))),i=1..9)},_C2)),x=0.85984.. .85986);**

This determines a value for .

`> `
**map(solve,subs(_C1=.85985,{seq(series_j(50,(1-(.5+.01*i))/9)=eval(subs(u=.5+.01*i,convert(op(2,dsol_ser_f),polynom))),i=1..9)}),_C2);**

`> `
**fsol:={_C1=.85985,_C2=convert(",`+`)/nops(")};**

Plugging the previous values into yields the following truncated asymptotic expansion of :

`> `
**subs(fsol,collect(subs(u=1-9*z,convert(series(op(2,dsol_ser_f),u,3),polynom)),[ln(1-9*z),z],normal));**

Once again, we compute a more precise expansion, but do not display it.

`> `
**j_num:=subs(fsol,collect(subs(u=1-9*z,convert(op(2,dsol_ser_f),polynom)),[ln(1-9*z),z],normal)):**

**Plotting **

We have three means to evaluate numerical values of : its Taylor series at 0, the differential equation it satisfies, the asymptotic expansion at .

The following plot is the numerically unstable result of following the graph of by using the differential equation. The method used is the fourth-order classical Runge-Kutta method.

`> `
**pl[RungeKutta]:=subs({LINE=POINT,COLOR(RGB,.9,.9,.2)=COLOR(RGB,0,0,0)},DEtools[DEplot](remove(type,double_inv_borel(power_j(3),j,z),equation),j(z),z=0.0001.. 1/9,[[j(.0000001)=1,D(j)(.0000001)=3]],j=0..5,stepsize=.0005)):**

The following plot is the result of the connection.

`> `
**pl[connection]:=plot(j_num,z=0..1/9,y=0..5,style=point,color=green):**

`> `
**plots[display](eval(subs(COLOUR=proc() COLOUR(RGB,1.,0,0) end,pl[Taylor])),pl[RungeKutta],pl[connection],title=`black=Runge-Kutta, red=Taylor, green=connection`);**

While the Taylor series converges rather slowly to the values of at the singularity, the numerical connection computed agrees with the curve obtained by following the differential equation. More specifically, the best approximation by the Taylor series is obtained for a truncation of order 100, whereas the expansion obtained by connection is of order 15 only.

We are now able to give an equivalent for the number of pairs of paths with same origin and end. Retaining only the singular part when tends to 0 of

`> `
**Order:=3: dsolve(d_eq,f(u),series);**

we get the equivalent
of
. An asymptotics for
follows by
*transfer of singularity*
: an asymptotic for the
-th coefficient a series is computed using the following table:

`> `
**matrix([[Function,`Asymptotic coefficient`,` `],[(1-z/rho)^(-alpha)*ln(1-z/rho)^beta,rho^n*n^(alpha-1)*ln(n)^beta/GAMMA(alpha),``(beta>=0 and alpha*` non integer`)],[(1-z/rho)^(-p)*ln(1-z/rho)^beta,rho^n*n^(p-1)*ln(n)^(beta-1)/GAMMA(alpha),``(beta>=0 and p*` integer`)],[1/ln(1-z/rho),rho^n/n/ln(n)^2,` `]]);**

In the case of , the result is

`> `
**asympt_p[3]:=-subs(fsol,_C2)*9^n/n;**

**Asymptotics for **

We now turn to our original goal: evaluating the coefficients of

.

It follows from the value of that the ordinary generating function behaves like where is close to and . Because of the positivity of the coefficients of , never vanishes before its singularity , so that also has no singularity before , as visualized on the next plot.

`> `
**plot(2/3-1/2/j_num,z=0..1/9,y=-.1..1,style=point,color=green);**

Correspondingly, we get the following asymptotics for , by transfer of singularity, using the table above:

`> `
**asympt_q[3]:=-subs(fsol,1/_C2)*9^n/n/ln(n)^2/2;**

**Numerical check of the asymptotics of **** and **

To get satisfactory results, we have to work for large .

`> `
**NMax:=250;**

We create another procedure to compute the coefficients of the series .

`> `
**proc_p:=rectoproc(diffeqtorec(double_inv_borel(power_j(3),j,z),j(z),u(n)),u(n),remember);**

We get a precise expansion for :

`> `
**series_p:=convert([seq(proc_p(i)*z^i,i=0..NMax)],`+`):**

from which follows a precise expansion for :

`> `
**series_q:=series(1/2+3*z/2-1/series_p/2,z,NMax+1):**

The asymptotic behaviour of is easily attained.

`> `
**for i from 50 to NMax by 50 do i=evalf(subs(n=i,asympt_p[3])/coeff(series_p,z,i)) od;**

The behaviour of , however, is not observed so obviously, due to a slow logarithmic convergence.

`> `
**for i from 50 to NMax by 50 do i=evalf(subs(n=i,asympt_q[3])/coeff(series_q,z,i)) od;**

Enlarging , of course, would lead to a limiting value of 1.

**Conclusion**

To conclude, we compute the table indicating the indicial polynomial at the dominant singularity and the possible asymptotic behaviours of there.

`> `
**behaviours:=proc(f,x,u) local L,i;
L:=sort([solve(f,x)]);
convert([seq(u^L[i]*ln(u)^nops(select(type,map(`+`,L[1..i-1],-L[i]),integer)),i=1..nops(L))],`+`)
end:**

`> `
**M[1]:=d,s,`indicial polynomial at s`,`possible behaviours for P[d](z)`:**

`> `
**for i from 2 to 7 do
M[i]:=i,1/i^2,
indicial_poly({algebraicsubs(double_inv_borel(power_j(i),j,z),j=(1-z)/i^2,j(z))},j,z);
B:=behaviours(M[i][3],eta,u);
M[i]:=M[i],B
od:**

`> `
**matrix([seq([M[i]],i=1..7)]);**

The function is a linear combination of the possible behaviours. Solving the connection would indicate which of them appear (with a non-zero coefficient) in . Next, transfer of singularity on the dominant behaviour would yield an asymptotic form for .