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

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

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:

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);

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)]);