Indicators of solvability for lattice models

A J Guttmann

Department of Mathematics & Statistics, The University of Melbourne Parkville, Vic. 3052, Australia

July 8, 1998

Abstract

We present a numerical method, based on exact series expansions, that distinguishes between lattice-based models both in combinatorics and statistical mechanics that are likely to be solvable in terms of standard functions of mathematical physics, and those that possess a natural boundary in a suitably defined complex plane. This latter class cannot therefore be algebraic, nor differentiably finite nor constructably differentiably algebraic. Known solutions in this latter case are all expressed as q-generalisations of standard functions.

1 Introduction

Some of the most famous results in mathematics involve a proof of the intrinsic unsolvability of certain problems. Some, such as ‘trisecting an angle’ are of long standing, while others, such as the lack of integer solutions to the equation x^{n} + y^{n} = z^{n} for n > 2 have only quite recently been acceptably proved [?]. In mathematical physics and combinatorics such results concerning the solvability or otherwise of problems are largely unknown. In this article we take a first step in addressing this absence by presenting and developing what is essentially a numerical method that provides, at worst, strong evidence that a problem has no solution within a large class of functions, including algebraic, differentiably finite (D-finite) [?, ?] and at least a sub-class [?] of differentially algebraic functions, called constructably differentiably algebraic (CDA). Since most of the special functions of mathematical physics — in terms of which most known solutions are given — are differentiably finite, this exclusion renders the problem unsolvable within this class. Throughout this article I will use the term D-unsolvable to mean that the problem has no solution within the class of D-finite functions as well as the sub-class of differentiably algebraic functions described above.

In fact, the exclusion is wider than this, as we show that the solutions possess a natural boundary on the unit circle in an appropriately defined complex plane. This excludes not only D-finite functions, but a number of others as well — though we have no simple way to describe this excluded class.

It may be worthwhile to recall the definitions of these classes of functions. A series f(z) is said to be differentiably finite if there exists an integer k and polynomials P_{0}(z), ⋯, P_{k}(z) with complex coefficients such that P_{0}(z) is not the null polynomial and

P_{0}(z)f(z) + P_{1}(z)f′(z) + ⋯ + P_{k}(z)f^{(k)}(z) = 0.

A series f(z) is said to be differentiably algebraic if there exists an integer k and a polynomial P in k + 1 variables with complex coefficients, such that

P (f(z), f′(z), ⋯, f^{(k)}(z)) = 0.

A series f(z) is said to be constructably differentiably algebraic if there exists both series f_{1}(z), f_{2}(z), ⋯, f_{k}(z) with f = f_{1}, and polynomials P_{1}, P_{2}, ⋯, P_{K} in k-variables, with complex coefficients, such that

f′_{1} = P_{1}(f_{1}, f_{2}, ⋯, f_{k}),  (1)
f′_{2} = P_{2}(f_{1}, f_{2}, ⋯, f_{k}),
⋯
f′_{k} = P_{k}(f_{1}, f_{2}, ⋯, f_{k}).

Differentiably algebraic functions in several variables are discussed in [?].

The method which we shall describe and which can, in favourable circumstances, be sharpened into a formal proof, has been applied to a wide variety of problems in both statistical mechanics and combinatorics. An underlying requirement is that the problem admits to a combinatorial formulation requiring the enumeration of graphs on a lattice. Typically, the solution of the problem will require the calculation of the graph generating function in terms of some parameter, such as perimeter, area, number of bonds or sites. A key first step is to anisotropise the generating function. For example, if counting graphs by the number of bonds on, say, an underlying square lattice, one distinguishes between horizontal and vertical bonds. In this way, one can construct a two-variable generating function, G(x, y) = ∑_{m,n} g_{m,n}x^{m}y^{n} where g_{m,n} denotes the number of graphs with m horizontal and n vertical bonds. Summing over one of the variables, we may write

G(x, y) = ∑_{m,n} g_{m,n}x^{m}y^{n} = ∑_{n} H_{n}(x)y^{n}  (2)

where H_{n}(x) is the generating function for the relevant graphs with n vertical bonds. It has been observed in all the problems so far studied, that the so called H functions are rational, with denominator zeros lying on the unit circle in the complex x plane.

In some cases one finds only a small finite number (typically one or two) of denominator zeros on the unit circle. Loosely speaking, this is the hallmark of a solvable problem. If, as is often observed, the denominator zeros become dense on the unit circle as n increases, so that in the limit a natural boundary is formed, then this is the hallmark of a D-unsolvable problem.

The significance of this observation is substantial. It is observed that, as n increases, the denominators of H_{n}(x) contain zeros given by steadily higher roots of unity. Hence the structure of the functions H_{n}(x) is that of a rational function whose poles all lie on the unit circle in the complex x-plane, such that the poles become dense on the unit circle as n gets large. This behaviour implies that the functions H_{n}(x) and hence G(x, y) as a function of x for y fixed, (a) have a natural boundary (b) are neither algebraic nor D-finite, nor CDA.

Of course, we are primarily interested in the solution of the isotropic case, when x = y, and it is clear that the anisotropic case can behave quite differently from the isotropic case. This is most easily seen by construction. Consider the function

f(x, y) = f_{1}(x, y) + (x − y)f_{2}(x, y),  (3)

where f_{1}(x, y) is D-finite and f_{2}(x, y) is not. Clearly, the function f(x, y) is not D-finite, while f(x, x) is D-finite. However, in all the cases we have studied where the solutions are known, the effect of anisotropisation does not change the analytic structure of the solution. Rather, it simply moves singularities around in the complex plane, at most causing the bifurcation of a real singularity into a complex pair. This can readily be seen from equation ??, given below, for the magnetisation of the Ising model. Replacing y by λx and varying λ merely causes the singularities to move smoothly, and indeed initially linearly, with λ in the complex plane. Further, for unsolved problems, numerical procedures indicate that a similar behaviour prevails. Nevertheless, this remains an observation, rather than an established fact, and, strictly speaking, should be established for each new problem.

If we now ask what functions do display the type of behaviour we have just observed - a build up of singularities on the unit circle in the complex plane, then the most obvious candidate that displays this behaviour is the qgeneralisation of the standard functions of mathematical physics. We have seen these in a number of solutions already, such as the hard hexagon model [?], certain interacting walk models [?] and some polygon models [?]. An explicit example is given immediately below.

That being said, not all problems with a small number of denominator zeros have been solved, while some D-unsolvable problems have been solved. In the former case however we believe that it is only a matter of time before a solution is found for these problems, while in the latter case the solutions have usually been expressed in terms of q-generalisations of the standard functions, which are of course not D-finite. As an example, the generating function for the number of parallelogram polygons given in terms of the area (q), horizontal semi-perimeter (x) and vertical semi-perimeter (y) is [?]

G(x, y, q) = y \frac{J_{1}}{J_{0}} where  (4)

J_{1}(x, y, q) = ∑_{n≥1} \frac{(−1)^{n−1}x^{n}q^{\binom{n+1}{2}}}{(q)_{n−1}(yq)_{n}} and  (5)

J_{0}(x, y, q) = ∑_{n≥0} \frac{(−1)^{n}x^{n}q^{\binom{n+1}{2}}}{(q)_{n}(yq)_{n}} .  (6)

In this case, it is clear that if we look at G(x, 1, q) in the complex q-plane with x held fixed, the solution possesses a natural boundary on the unit circle.

We suggest that this method is a particularly useful first step in the study of such problems. One anisotropises, generates enough terms in the generating function to be able to construct the first few H functions, then studies the denominator pattern. If it appears that the zeros are becoming dense on the unit circle, one has good reason to believe that the problem is D-unsolvable. If on the other hand there are only one or two zeros, one is in an excellent position to seek the solution in terms of the standard functions of mathematical physics — loosely defined as those described in [?]. In some cases one may be able to prove that the observed denominator pattern persists. In that case, one has proved the observed results.

The construction of the H functions deserves some explanation. At very low order this can often be done exactly, by combinatorial arguments based on the allowed graphs. Beyond this, our method is to generate the coefficients in the expansion, assume it is rational, then by essentially constructing the Padé approximant one can conjecture the solution. Typically, one might generate 50-100 terms in the expansion and find a rational function with numerator and denominator of perhaps degree 5 or 10. Thus the first 10 or 20 terms of the series are used to identify the rational function, the remainder are used to confirm it. Hence while this is not a derivation that proves that the function is rational, the chance of it being irrational is extraordinarily small.

It should be said explicitly that this technique is computationally demanding. That is to say, the generation of sufficient terms in the generating function is usually quite difficult. Only with improved algorithms — most notably the combination of the finite lattice method [?, ?] with a transfer matrix formulation — and computers with large physical memory that are needed for the efficient implementation of such algorithms, has it been possible to obtain expansions of the required length in a reasonable time. The technique is still far from routine, with each problem requiring a significant calculational effort.

An additional, and exceptionally valuable feature of the method comes when the numerical work, described above, is combined with certain functional relations that the anisotropised generating functions must satisfy. In the language of statistical mechanics, these key functional relations are called inversion relations and imply a connection between the generating function and its analytic continuation, usually involving the reciprocal of one or more of the expansion variable(s). As we show below, the existence of these inversion relations, coupled with any obvious symmetries (usually a symmetry with respect to the interchange of x and y), coupled with the observed behaviour of the H functions — described above — can yield an implicit solution to the underlying problem with no further calculation. An example of this is the solution [?] of the zero-field free energy of the two-dimensional Ising model.

In the remainder of this article, we describe the method in considerable detail in a few cases, then go on to apply it to a range of problems in statistical mechanics and combinatorics. We also take the first steps in extending the inversion relation idea from its natural home in statistical mechanics to the arena of combinatorics — where it sits less naturally due to the absence of an underlying Hamiltonian, the symmetries of which give rise to the inversion relation. It is comforting to discover that, without exception, the unsolved problems we discuss are all found to be D-unsolvable.

Other important aspects of the method, such as the connection of these ideas with concepts of integrability, and with the existence of a Yang-Baxter equation are not explored here.

2 The Ising model free energy and magnetisation

The one-dimensional Ising model consists of a chain of N spins, each of which may point up or down, denoted μ_{i} = ±1, i = 1, ⋯ N. Each spin interacts only with nearest-neighbour spins with interaction strength J and with an external magnetic field, with interaction strength H. This model is described [?] by the Hamiltonian

H = −J ∑_{<i,j>} μ_{i}μ_{j} − H ∑_{i=1}^{N} μ_{i}  (7)

where the first sum is over nearest-neighbour pairs. Imposing cyclic boundary conditions, so that μ_{N+1} = μ_{1} allows us to write the first sum explicitly as ∑_{i=1}^{N} μ_{i}μ_{i+1}.

The partition function, in the thermodynamic limit, is defined by

lim_{N→∞} Z(N, T)^{\frac{1}{N}} ,

where

Z(N, T) = ∑_{μ_{1}=±1} ∑_{μ_{1}=±2} ⋯ ∑_{μ_{N} =±1} exp(−βH)

and β = 1/k_{B}T where k_{B} is Boltzmann’s constant.

This multiple sum can be expressed as an iterated matrix product [?] and the problem then reduces to finding the eigenvalues of a 2 × 2 matrix.

The result is

Z(K, H) = exp(K) cosh H + \sqrt{exp(2K) sinh^{2} H + exp(−2K)} (8)

where K = βJ. It can readily be verified that the partition function satisfies a so called inversion relation

Z(K, H)Z(K + \frac{iπ}{2} , −H) = 2i sinh(2K),  (9)

which connects the partition function and its analytic continuation.

In two dimensions the problem is substantially more difficult [?]. If we take a lattice of M rows and N columns, then the first term in the Hamiltonian now becomes a double sum

J_{1} ∑_{i=1}^{M−1} ∑_{j=1}^{N} μ_{i,j}μ_{i+1,j} + J_{2} ∑_{i=1}^{M} ∑_{j=1}^{N} μ_{i,j}μ_{i,j+1},  (10)

where cylindrical boundary conditions have been imposed, so that μ_{i,N+1} = μ_{i,1}, i = 1, 2, ⋯, M.

The calculation of the partition function now involves the diagonalization of a 2^{M} × 2^{M} matrix in the limit as M → ∞, a calculation which has only been carried out [?], [?] in the case of zero magnetic field (H = 0). In the limit of an infinitely large lattice (lim_{M,N} → ∞) one finds

log Z(K_{1}, K_{2}) = log 2 + \frac{1}{2π^{2}} ∫ ∫_{0}^{π} log f(θ_{1}, θ_{2})dθ_{1}dθ_{2},

where K_{1} = βJ_{1}, K_{2} = βJ_{2} and

f(θ_{1}, θ_{2}) = cosh 2K_{1} cosh 2K_{2} − sinh 2K_{1} cos θ_{1} − sinh 2K_{2} cos θ_{2}

which can be expressed in terms of complete elliptic integrals. This is one of the most famous calculations of 20^{th} century statistical mechanics.

It is convenient to define the reduced partition function by

Λ(t_{1}, t_{2}) = Z(K_{1}, K_{2})/2 cosh K_{1} cosh K_{2},

where t_{1,2} = tanh K_{1,2}. The reduced partition function then [?] satisfies the inversion relation

ln Λ(t_{1}, t_{2}) + ln Λ(1/t_{1}, −t_{2}) = ln(1 − t_{2}^{2}),

where again the second term in the sum is an analytic continuation of the first. Writing the reduced partition function

ln Λ(t_{1}, t_{2}) = ∑_{n,m} a_{n,m}t_{1}^{2m}t_{2}^{2n} = ∑_{n} R_{n}(t_{1}^{2})t_{2}^{2n}

Baxter [?] points out that Onsager’s solution can be used to show that

R_{n}(t_{1}^{2}) = P_{2n−1}(t_{1}^{2})/(1 − t_{1}^{2})^{2n−1},

where P_{2n−1}(t^{2}) is a polynomial in t^{2} of degree 2n−1. That is, the functions R_{n} are rational, with numerator and denominator of equal degree, and with the denominator having only one pole of degree 2n − 1 in the complex t_{1}^{2} plane, at t_{1}^{2} = 1.

The significance of this observation is that when it is coupled with the above inversion relation and the obvious symmetry

Λ(t_{1}, t_{2}) = Λ(t_{2}, t_{1}),

it is sufficient to determine, order by order, the numerator polynomials P_{n}. That is to say, the complete Onsager solution is implicitly determined by these two functional equations and the simple form of the denominator (and some analyticity assumptions [?]).

To be more precise, as the numerator is a polynomial of degree 2n − 1 there are 2n unknowns in the numerator. The symmetry relation reduces this to n unknowns, and the inversion relation allows us to find the n unknowns.

A similar result is seen if we consider the spontaneous magnetisation of the anisotropic square lattice Ising model. Denoting x = exp(−4J_{1}/k_{B}T), y = exp(−4J_{2}/k_{B}T), the magnetisation is [?]

M(x, y) = [1 − \frac{16xy}{(1 − x)^{2}(1 − y)^{2}} ]^{\frac{1}{8}} .  (11)

This clearly satisfies the symmetry relation M(x, y) = M(y, x), and it also satisfies the inversion relation M(x, y)−M(x, 1/y) = 0, as can be seen upon inspection. Writing the magnetisation as a partially re-summed series, so that

M(x, y) = 1 − ∑_{n} H_{n}(y)x^{n},

we observe that the functions H_{n} are rational functions of the form

H_{n}(y) = \frac{yP_{n}(y)}{(1 − y)^{2n}} ,

for |y| < 1 where P_{n}(y) is a polynomial of degree 2n − 2. Each such function can be analytically continued to |y| > 1, and substitution into the inversion relation allows one to verify that it is satisfied, as far as one cares to push the expansion.

The first few numerators are

−2y , −2y(2 + 3y + 2y^{2}) , −2y(3 + 16y + 32y^{2} + 16y^{3} + 3y^{4}), for n = 1, 2, 3.

As observed for the free energy, the symmetry relation, inversion relation and functional form of the H functions are sufficient to determine the solution. Some other models were similarly solved by Stroganov [?].

Note that two independent features were necessary for this method of solution. The existence of the inversion and symmetry relation is one feature, and the particularly simple form of the functions R_{n}, and in particular their denominator structure is the other. These two examples led us to try and generalise this approach to other solved and unsolved problems, in order to obtain solutions, or at least additional information.

For statistical mechanical systems, the existence of an inversion relation follows from the underlying symmetries of the Hamiltonian [?, ?, ?, ?, ?]. As a consequence, a number of unsolved problems, such as the susceptibility of the two-dimensional Ising model, the free energy of the three-dimensional Ising model, and various thermodynamic properties of the q-state Potts model [?] all possess known inversion relations. For example, for the free energy of the three-dimensional Ising model [?], if the reduced partition function is defined by Λ(t_{1}, t_{2}, t_{3}) = Z(K_{1}, K_{2}, K_{3})/2 cosh K_{1} cosh K_{2} cosh K_{3}, where t_{1,2,3} = tanh K_{1,2,3}, the reduced partition function then satisfies the inversion relation

ln Λ(t_{1}, t_{2}, t_{3}) + ln Λ(1/t_{1}, −t_{2}, −t_{3}) = ln(1 − t_{2}^{2}) + ln(1 − t_{3}^{2}).

For lattice based combinatorial structures there is in general no analogue of a Hamiltonian, and so notions of the symmetry group of the Hamiltonian are not relevant. There is thus no obvious concept of an inversion relation. However it is possible to determine analogous functional relations for some combinatorial problems, though to date these have largely been derived “experimentally”, as we show below.

As well as the inversion relation (and symmetry relation), the general form of the rational functions (denoted R_{n} in the case of the Ising partition function above) needs to be known. Unless the problem has already been fully solved, this will not usually be a priori known. In fact it too will be determined “experimentally” and, in favourable cases, subsequently proved.

These observations lead to the following proposed approach to the study of statistical mechanical systems for which inversion relations are known. We derive the (anisotropic) series expansion of the quantity of interest, partially re-sum as above, and study the analytic properties of the functions which are the coefficients of the re-summed series. As mentioned above, the derivation of the series is usually a demanding computational exercise, for which efficient algorithms need to be designed. Otherwise there is simply insufficient data for the above approach to be pursued. It is the development of such algorithms and the availability of cheap, fast computing that has made this approach possible.

3 The Ising model susceptibility.

As our first example of this proposed approach to the study of unsolved problems, we consider the susceptibility of the two-dimensional Ising model, which is one of the most extensively studied [?], yet still unsolved, problems in statistical mechanics. For the square lattice version of this model, the relevant inversion relation [?] is χ(t_{1}, t_{2}) + χ(1/t_{1}, −t_{2}) = 0, and the symmetry relation χ(t_{1}, t_{2}) = χ(t_{2}, t_{1}) also holds. Partially 2r allows us to write the susceptibility as

χ(t_{1}, t_{2}) = ∑_{m,n≥0} c_{m,n}t_{1}^{m}t_{2}^{n} = ∑_{n≥0} H_{n}(t_{1})t_{2}^{n}.

This approach was first taken in [?], in which H_{1}, H_{2}, H_{3}, and H_{5} were found. They are

H_{0}(t) = (1 + t)/(1 − t),
H_{1}(t) = 2(1 + t)^{2}/(1 − t)^{2},
H_{2}(t) = 2(1 + 6t + 8t^{2} + 6t^{3} + t^{4})/(1 − t)^{3}(1 + t),
H_{3}(t) = 2(1 + 8t + 10t^{2} + 8t^{3} + t^{4})/(1 − t)^{4} and
H_{5}(t) = 2(1 + 16(t + t^{7}) + 64(t^{2} + t^{6}) + 144(t^{3} + t^{5}) + 166t^{4} + t^{8})/(1 − t)^{6}(1 + t)^{2}.

The calculation of these functions is computationally demanding, being of exponential complexity. Over the last twenty years, Enting [?] has developed an alternative method, known as the finite lattice method, which while still of exponential complexity, is nevertheless exponentially faster than preexisting methods based on direct enumeration. Based on this method, we [?] have obtained the first 14 of these rational functions.

Even from the 5 H functions given above, it is clear that the situation is not as straightforward as for the partition function or magnetization. For the next two, we [?] find

H_{4}(t) = 2(1 + t^{10} + 15(t + t^{9}) + 71(t^{2} + t^{8}) + 192(t^{3} + t^{7}) + 326(t^{4} + t^{6}) + 388t^{5})
/(1 − t^{3})(1 − t)^{4}(1 + t)^{3},
H_{6}(t) = 2(1 + t^{18} + 28(t + t^{17}) + 220(t^{2} + t^{16}) + 1149(t^{3} + t^{15}) + 4081(t^{4} + t^{14})

+10768(t^{5} + t^{13}) + 22083(t^{6} + t^{12}) + 36283(t^{7} + t^{11}) + 48543(t^{8} + t^{10})

+53446t^{9})/(1 − t^{3})^{3}(1 − t)^{4}(1 + t)^{5}.

For all n ≤ 14 (including the others not shown), the numerator polynomial is found to be symmetric, unimodal and with positive coefficients. The denominator polynomial has zeros lying on the unit circle, at t = 1 for all n, at t = −1 for n = 2 and n ≥ 4, and we observe that for n = 4 and n ≥ 6 there are zeros at t^{3} = 1 and for n = 12 and n ≥ 14 there are zeros at t^{5} = 1. The numerator and denominator are of equal degree, notably 1, 2, 4, 4, 10, 8, 18, 20, 26, 28, 34, 36, 48, 44, 62 ⋯ for n = 0, 1, 2, ⋯, 14 respectively.

The degree of the polynomials is increasing so rapidly that even if we could predict the denominator for all n, the constraints imposed by the inversion relation and the symmetry relation are insufficient to implicitly yield the solution, unlike the case of the free energy and magnetisation.

Nevertheless, we can obtain useful analytic information about the structure of the solution. We use the notation of [?], in which the susceptibility of the Ising model is expressed as an expansion in terms of so called 2k + 1 particle excitations,

χ(t_{1}, t_{2}) = ∑_{k} χ^{2k+1}(t_{1}, t_{2}).

We note in passing that this expansion applies to the high-temperature susceptibility. For the low temperature susceptibility the corresponding expansion involves 2k particle excitations. The notion of particle excitations is the language of a field-theoretic expansion of the Ising model, an explanation of which would take us unnecessarily far afield. It suffices to say that such an expansion exists, and refer the interested reader to [?] for details.

Syozi and Naya [?] appear to have been the first to calculate χ^{1}, even though their calculation preceded the particle excitation formulation of [?]. ¿From [?] we find

χ^{1} = \frac{(1 − t_{1}^{2})(1 − t_{2}^{2})}{(1 − t_{1} − t_{2} − t_{1}t_{2})^{2}} (1 − 16\frac{t_{1}^{2}t_{2}^{2}}{(1 − t_{1}^{2})^{2}(1 − t_{2}^{2})^{2}} )^{\frac{1}{4}} = ∑_{n} H_{n}^{(1)}(t_{1})t_{2}^{n},

and H_{n}^{(1)}(t_{1}) = H_{n}(t_{1}) for n = 1, 2, 3, 5 while

H_{4}^{(1)}(t) = 2(1+t^{8}+14(t+t^{7})+56(t^{2}+t^{6})+122(t^{3}+t^{5})+146t^{4})/(1−t)^{5}(1+t)^{3}.

It is straightforward to show that, for all n, the numerators are symmetric, unimodal polynomials (with positive coefficients). Further, for n even, the denominator is (1 − t)^{n+1}(1 + t)^{n−1}, and for n odd, the denominator is (1 − t)^{n+1}(1 + t)^{n−3}. In both cases negative subscripts are to be replaced by zero. The structure of the numerator and denominator, taken together, imply that the symmetry and inversion relations that hold for χ also hold for χ^{1}. For n even, the numerator and denominator polynomials are of degree 2n, hence there are 2n + 1 coefficients to be determined. Symmetry reduces this to n + 1, while the inversion relation determines n coefficients, leaving 1 unknown. This can be determined by the observation that the residue at t = −1 of H_{n}^{(1)}(t), for n even, is

−\frac{(2n − 5)!!!!}{2(n/2)!}

where n!!!! = n(n − 4)(n − 8) ⋯, terminating at the smallest integer greater than 0. For n odd, the numerator and denominator polynomials are of degree 2n − 2, hence there are 2n − 1 coefficients to be determined. Symmetry reduces this to n − 1, and the inversion relation determines all of these.

For χ^{3}, χ^{5}, etc. no closed form expression is yet known - though they can [?] be expressed as hyper-elliptic integrals, and at least the first few are [?] differentiably finite. (They probably all are but this hasn’t been proved.) However our numerical studies clearly imply (but do not prove) that a similar, but more complex structure prevails in these cases.

The principal features we observe are that χ^{2k+1} can be similarly expanded in terms of rational functions, as shown explicitly for χ^{1} above, with numerators and denominators of equal degree. Furthermore, the numerators are observed to be unimodal, symmetric and with all coefficients positive, from which follows that the symmetry and inversion relations apply not only to χ, but to each term χ^{2k+1} in its expansion — at least as far as we have proceeded.

Further, we find that the denominator of the rational coefficients which occur in the expansion of χ^{2k+1} have, in addition to the factors in χ^{1} given above, systematic occurrences of powers of the terms (1−t^{3}), (1−t^{5}), ⋯, (1− t^{2k+1}), which can be predicted [?]. ¿From the results for χ and χ^{1} given above, it can be seen that the first contribution of χ^{3} to χ occurs in H_{4}, (as evidenced by the occurrence of the term (1 − t^{3}) in the denominator). Similarly, we find that the first contribution of χ^{5} to χ occurs in H_{12}. Hence it appears that the first occurrence of the factor 1−t^{2k+1} in the denominator coincides with the first contribution of a 2k + 1 particle excitation.

It follows that, as n increases, the denominators of H_{n}(t) contain zeros given by the (2k + 1)^{th} roots of unity. And as n increases, so does k. Hence the structure of the functions H_{n}(t) is that of a rational function whose poles all lie on the unit circle in the complex t-plane, such that the poles become dense on the unit circle as n gets large. This behaviour implies (unless miraculous cancellation of almost all poles suddenly starts to occur at high order) that the functions H_{n}(t) and hence χ(t_{1}, t_{2}) as a function of t_{1} for t_{2} fixed, (a) have a natural boundary (b) are neither algebraic nor D-finite, nor CDA.

Leaving these considerations aside for the moment, the significance of these observations of the Ising model is that the observed behaviour suggests a new and powerful tool to investigate the analytic structure of a wide variety of problems. By generalizing to the anisotropic model, and studying the distribution of zeros of the denominators in the H_{n}(t) functions and their analogues, we can distinguish between those that are likely to be solvable in terms of standard functions, and those that are not. In the former case there is a finite number (usually one or two) of singularities on the unit circle, while in the latter case there is, in the limit of large n, an infinite number, corresponding to a natural boundary. Numerically, this is signified by an increasing number of singularities in the denominator of H_{n}(t) as n increases. In favourable cases one can predict the behaviour of the denominator, and thus prove that the number of denominator zeros grow indefinitely. The magnetisation and susceptibility of the two-dimensional zero field Ising model, discussed above, are examples of the two types of behaviour.

We do not claim that in the former case the solution is D-finite, though all the models we have studied that display this behaviour are D-finite. Indeed, it is easy to construct an example of non D-finite functions that display this behaviour. For example,

f(x, y) = exp(x(exp \frac{y}{1 − x}−1)) = 1+\frac{xy}{1 − x}+\frac{x(1 + x)y^{2}}{2(1 − x)^{2}} +\frac{x(1 + 3x + x^{2})y^{3}}{6(1 − x)^{3}} +· · · .

(12)

Another example, corresponding to a model of interest, rather than just a contrived function as above, is that of three-dimensional convex polygons [?, ?] enumerated by perimeter. The generating function is not D-finite, yet the H-functions appear to have only a single denominator zero (though this has not been tested at high order).

A related but distinct observation is that the existence of inversion relations, coupled with construction of the H-functions, provides an alternative method of solution in some cases. We show, in the next section, how this concept can be applied to certain combinatorial problems too.

4 Staircase polygons

The enumeration of staircase polygons by perimeter is one of the simpler combinatorial exercises, but is nevertheless useful pedagogically, as so many distinct methods can be demonstrated in its solution. To this long list we add the experimental approach of studying the early terms of the two variable series expansion of the perimeter generating function and observing a functional relation equivalent to the inversion relation discussed above for certain statistical mechanical systems.

If we write the perimeter generating function as

P (x, y) = ∑_{n,m} p_{n,m}x^{2n}y^{2m} = ∑_{m} H_{m}(x^{2})y^{2m}

then H_{m}(x^{2}) is the generating function for staircase polygons with 2m vertical bonds.

¿From observation of the early terms, it is clear that

H_{m}(x^{2}) = x^{2}S_{m}(x^{2})/(1 − x^{2})^{2m−1}

for m > 1, where S_{m}(x^{2}) is a symmetric, unimodal polynomial with nonnegative coefficients, of degree (m − 2). This observed symmetry can be expressed formally as

x^{2m}H_{m}(x^{2}) + x^{2}H_{m}(1/x^{2}) = 0, n > 1.

This in turn translates into the functional relation

P (x, y) + x^{2}P (1/x, y/x) = −y^{2}.

There is also an obvious symmetry relation P (x, y) = P (y, x), and these observations are sufficient to implicitly solve the problem by calculating the H functions order by order in polynomial time.

Of course, this must rank as one of the least impressive ways of solving this fairly simple model. However the purpose of this example is twofold. Firstly to show that this essentially experimental method can be applied to combinatorial structures in order to discover an inversion relation. Secondly, to show that once one has such an inversion relation, then this, coupled with symmetry and the structure of the H functions, (plus certain analyticity assumptions) provides an alternative method for obtaining a solution (albeit experimentally). Once one has such a conjectured solution, it is a comparatively easy task to prove that it is correct.

Numerous other polygon problems can also be tackled similarly [?].

5 Three-choice polygons

The problem of three-choice polygons [?] is an intriguing one, as we know everything about this model except a closed form solution. We have a polynomial time algorithm to generate the coefficients in its series expansion — which is tantamount to a solution — and have made an analysis of its asymptotic behaviour.

We have recently anisotropised the model [?, ?] in order to see whether the ideas developed here give insight into the solution.

Let P_{3}(x, y) = ∑_{m,n} a_{m,n}x^{m}y^{n} be the perimeter generating function, where a_{m,n} gives the number of 3-choice polygons, distinct up to translation, with 2m horizontal bonds and 2n vertical bonds. Then

P_{3}(x, y) = ∑_{n} H_{n}(y)x^{n},

where

H_{n}(y) = P_{n}(y)/Q_{n}(y)

is a rational function of y. The degree of the numerator polynomial increases like 3n while the denominators are observed to be

Q_{n}(y) = (1 − y)^{2n−1}(1 + y)^{2n−7}, n even,
= (1 − y)^{2n−1}(1 + y)^{2n−8}, n odd,

where there are no terms in (1 + y) for n < 5. It is not difficult to construct a combinatorial argument, based on the way the polygon can “grow” that is consistent with this behaviour. This argument has recently been sharpened to a proof [?]. It has also been proved [?] that the solution is D-finite, and cannot be algebraic as the asymptotic behaviour of the number of coefficients [?] includes a logarithmic term.

An inversion relation for this model can be found experimentally [?], and the solution possesses (x, y) symmetry. Nevertheless, because the numerator grows like 3n we do not have enough constraints to implicitly solve the model. What is needed is some additional constraint on the behaviour of the coefficients, the discovery of which has so far eluded us. We nevertheless consider this a promising approach, which has already revealed valuable analytic information about the solution.

6 Hexagonal directed animals

A directed site animal A on a lattice is defined to be a set of vertices such that all vertices p ∈ A are either the (unique) origin vertex or may be reached from the origin by a connected path, containing bonds only in the allowed directions, through sites of A.

In [?, ?] it was found that the number of such animals of perimeter n grew asymptotically like μ^{n}/\sqrt{n}, where μ = 4 for the triangular lattice, μ = 3 for the square lattice. Furthermore, the generating function was given by the solution of a simple algebraic equation. For the hexagonal lattice however we [?] found μ = 2.025131 ± 0.000005, and were unable to solve for the generating function.

In order to gain more insight into this seemingly anomalous situation, the model was anisotropised [?]. Let A_{h}(x, y) = ∑_{m,n} a_{m,n}x^{m}y^{n} be the site generating function, where a_{m,n} gives the number of hexagonal lattice site animals, with m sites supported [?] one particular way and n − m other sites. Then

A_{h}(x, y) = ∑_{n} H_{n}(y)x^{n},

where

H_{n}(y) = P_{n}(y)/Q_{n}(y)

is observed to be a rational function of y.

However the denominator pattern, while regular, contains terms of the form (1 − y^{k}) where k is an increasing function of n. The degree of the numerator also increases faster than linearly. Thus this model displays the same qualitative behaviour as the susceptibility of the Ising model, discussed above. Hence similar conclusions such as the existence of a natural boundary in the appropriate complex plane, and that the solution is not D-finite may be drawn.

This is then consistent with the seemingly anomalous value of the constant μ.

7 Self avoiding walks and polygons

The problem of square lattice self avoiding walks (SAW) and self avoiding polygons (SAP) are much studied problems, equally widely known for their mathematical interest and their intractability. See for example [?, ?].

A study of anisotropic square lattice SAW has been reported in [?]. Writing the SAW generating function C(x) in the by now familiar form as

C(x, y) = ∑_{m,n≥0} c_{m,n}x^{m}y^{n} = ∑_{n≥0} H_{n}(x)y^{n},  (13)

we found the first 11 H functions, H_{0}(x), ⋯, H_{10}(x).

The first few are:

H_{0}(x) = (1 + x)/(1 − x),
H_{1}(x) = 2(1 + x)^{2}/(1 − x)^{2},
H_{2}(x) = 2(1 + 7x + 14x^{2} + 16x^{3} + 9x^{4} + 3x^{5})/(1 − x)^{3}(1 + x)^{2} and
H_{3}(x) = 2(1 + 10x + 29x^{2} + 44x^{3} + 41x^{4} + 22x^{5} + 7x^{6})/(1 − x)^{4}(1 + x)^{2}.

The first occurrence of the term (1 − x^{3}) appears in H_{5}(x) and the term (1 + x^{2}) first appears in H_{7}. Higher order roots of ±1 then systematically occur as n increases. The denominator pattern appears to be predictable, though we have not been able to prove this. The degree of the numerator is equal to the degree of the denominator in all cases observed.

Thus we see again the, by now, characteristic hallmark of a D-unsolvable problem. Similar behaviour is observed for SAP. Writing the SAP generating function P (x) as

P (x, y) = ∑_{m,n≥1} p_{m,n}x^{2m}y^{2n} = ∑_{n≥1} R_{n}(x)y^{2n},

where p_{m,n} is the number of square lattice polygons, equivalent up to a translation, with 2n horizontal steps and 2m vertical steps. We [?] found the first 9 R functions, R_{1}(x), ⋯, R_{9}(x), and these were found to behave in a manner characteristic of D-unsolvable problems — that is, the zeros appear to build up on the unit circle. They first few are:

R_{1}(x) = x/(1 − x),
R_{2}(x) = x(1 + x)^{2}/(1 − x)^{3},
R_{3}(x) = x(1 + 8x + 17x^{2} + 12x^{3} + 3x^{4})/(1 − x)^{5},
R_{4}(x) = x(1 + 18x + 98x^{2} + 204x^{3} + 178x^{4} + 70x^{5} + 11x^{6})/(1 − x)^{7},
R_{5}(x) = xP_{9}(x)/(1 − x)^{9}(1 + x)^{2},
R_{6}(x) = xP_{15}(x)/(1 − x)^{11}(1 + x)^{4},
R_{7}(x) = xP_{20}(x)/(1 − x)^{13}(1 + x)^{6}(1 + x^{2} + x^{4}).

In the above equations, P_{k}(x) denotes a polynomial of degree k. As for 3-choice polygons, the occurrence of new terms in the denominator, corresponding to higher roots of unity, can be identified with the first occurrence of specific graphs. In this way [?] the denominator pattern can be predicted, though rather more work is required to refine this observation into a proof.

A similar study of hexagonal lattice polygons [?] leads to similar conclusions. Furthermore, we observed that the denominators of the H functions for the square and hexagonal lattices are simply related.

8 The 8-vertex model

Very recently, as a test of the idea that “solvable” models should, when anisotropised, have H functions with only one or two denominator zeros, Tsukahara and Inami [?] studied the 8-vertex model - which is one of the most difficult statistical mechanics models that has been exactly solved [?]. It can be described as two inter-penetrating planar Ising models, coupled by a four-spin coupling, with two spins in each of the sub-lattices. Let the coupling in one sub-lattice be L, that in the other be K, and the four-spin coupling be M. Then the usual high-temperature expansion variables are t_{1} = tanh K, t_{2} = tanh L, t_{3} = tanh M. New high temperature variables may be defined as follows:

z_{1} = \frac{t_{1} + t_{2}t_{3}}{1 + t_{1}t_{2}t_{3}} ,  (14)
z_{2} = \frac{t_{2} + t_{1}t_{3}}{1 + t_{1}t_{2}t_{3}} ,  (15)
z_{3} = \frac{t_{3} + t_{1}t_{2}}{t_{1} + t_{2}t_{3}} .  (16)

then it has recently been shown [?] that the logarithm of the reduced partition function per face

log Λ(z_{1}, z_{2}, z_{3}) = ∑_{l,m,n} a_{l,m,n}z_{1}^{2l}z_{2}^{2m}z_{3}^{2n}

satisfies

log Λ(z_{1}, z_{2}, z_{3}) + log Λ(\frac{1 − z_{2}^{2}}{z_{1}(1 − z_{3}^{2})}, −z_{2}, −z_{3}) = log(1 − z_{2}^{2}).  (17)

A partial re-summation allows the reduced partition function to be written as

log Λ(z_{1}, z_{2}, z_{3}) = ∑_{m,n} R_{m,n}(z_{1}^{2})z_{2}^{2m}z_{3}^{2n}.  (18)

After a complicated graphical calculation [?], it was found that

R_{1,0}(z^{2}) = z^{2}/(1 − z^{2}),  (19)
R_{1,1}(z^{2}) = 2z^{4}/(1 − z^{2})^{3},  (20)
R_{2,0}(z^{2}) = z^{2}(2 − 5z^{2} + z^{4})/(1 − z^{2})  (21)
R_{1,2}(z^{2}) = 3z^{6}(1 + z^{2})/(1 − z^{2})^{5}.  (22)

It is then argued [?] that the general form of the coefficients is

R_{m,n}(z^{2}) = P_{m,n}(z^{2})/(1 − z^{2})^{2m+2n−1}.

This behaviour then accords with the expected behaviour of solvable models. That is to say, only a finite number — in this case 1 — of denominator singularities.

9 The three-dimensional Ising model

These ideas are also applicable to three-dimensional models, such as the three-dimensional Ising model. For this model there are no exact results known. However inversion relations can still be proved, (indeed, the appropriate relation in the case of the free energy is given in Section ??). Similarly, the susceptibility of the model on the simple cubic lattice, anisotropic in all three directions, satisfies the inversion relation

χ(t_{1}, t_{2}, t_{3}) + χ(1/t_{1}, −t_{2}, −t_{3}) = 0.  (23)

Here t_{1} = tanh(J_{1}/k_{B}T), t_{2} = tanh(J_{2}/k_{B}T), and t_{3} = tanh(J_{3}/k_{B}T). Partially re-summing allows us to write the susceptibility as

χ(t_{1}, t_{2}, t_{3}) = ∑_{l,m,n≥0} c_{l,m,n}t_{1}^{l}t_{2}^{m}t_{3}^{n} = ∑_{m,n≥0} H_{m,n}(t_{1})t_{2}^{m}t_{3}^{n}.  (24)

This approach was first taken in [?], in which H_{1,1}, H_{2,1}, H_{3,1}, and H_{2,2} were studied, though only the first three were identified. They were found to be

H_{1,1}(t) = 8(1 + t)^{3}/(1 − t)^{3},  (25)
H_{2,1}(t) = 16(1 + 5t^{2} + 7t^{2} + 5t^{3} + t^{4})/(1 − t)^{4},  (26)
H_{3,1}(t) = 8(3 + 23t + 46t^{2} + 46t^{3} + 23t^{4} + 3t^{5})/(1 − t)^{5}. (27)

These display the simple behaviour also observed for the first few Hfunctions in the case of the two-dimensional model. That is to say, the denominator has only a single zero. However the next term, which we have managed to identify from the raw data given in [?], already shows the occurrence of a cube root of unity in the denominator. It is

H_{2,2}(t) = 16(3 + 3t^{10} + 34(t + t^{9}) + 143(t^{2} + t^{8}) + 373(t^{3} + t^{7})
+623(t^{4} + t^{6}) + 745t^{5})/(1 − t^{3})(1 − t)^{4}(1 + t)^{3}.  (28)

Indeed, in structure it is very similar to H_{4}(t) for the two-dimensional case. We have not gone further, but it seems very likely that higher order H functions will be rational with denominators corresponding to higher roots of unity.

10 Percolation and directed percolation

Two widely studied but unsolved problems are that of ordinary and directed percolation in dimension two and higher. Our preliminary studies [?] show that these two models display the characteristic behaviour of D-unsolvable models — a build up of higher roots of unity in the denominators of the H functions when the models are anisotropised. Much further work needs to be done for these two models.

11 Conclusion

At worst we have developed a powerful numerical technique capable of indicating whether a problem is likely to be readily D-solvable or not. The hallmark of unsolvability, which is the build up of zeros around the unit circle in the complex x-plane in the H(x) functions of the anisotropised models, can, in favourable cases, be refined into a proof.

In the most favourable cases, where in addition an inversion relation can be obtained — as in the case of the free energy of the two-dimensional Ising model — or numerically inferred, as in the case of staircase polygons, and the H functions are sufficiently simple, with only a pole at 1 on the unit circle, an exact solution can be implicitly obtained.

Most tantalisingly, the prospect of solving hitherto unsolved problems by predicting the numerator and denominator of the H functions by a combination of combinatorial and symmetry based arguments remains open.

Clearly, much further work remains to be done in classifying precisely what class of functions is excluded by certain observed behaviour, and in developing methods to solve problems which are identifiable as D-unsolvable.

12 Acknowledgments

I would like to thank my colleagues Mireille Bousquet-Mélou, Richard Brak, Will Orrick, Paul Pearce and Aleks Owczarek for valuable comments on the manuscript, and Mireille Bousquet-Mélou and Jean-Marie Maillard for many illuminating discussions. Much of the underlying enumerative work was carried out in collaboration with Ian Enting. This work was supported by the Australian Research Council.

References

[1] M. Abramowitz and I. A. Stegun (Editors), Handbook of Mathematical Functions, (Dover, New York) (1972).

[2] R. J. Baxter, in Fundamental Problems in Statistical Mechanics, 5, E. G. D. Cohen ed. (North Holl. Amsterdam) (1981).

[3] R. J. Baxter, Hard hexagons: exact solution, J. Phys. A. 13, L61-70, (1980).

[4] R. J. Baxter, The Inversion relation method for some two-dimensional exactly solved models in lattice statistics, J. Stat. Phys, 28, 1-41 (1982).

[5] M. Bousquet-Mélou, A method for the enumeration of various classes of column-convex polygons, Disc. Math. 154, 1-25, (1996).

[6] M. Bousquet-Mélou, (Private communication).

[7] M. Bousquet-Mélou and X. Viennot, Empilements de segments et qénumération de polyominos convexes dirigés. J. Combin. Theory Ser. A 60), (1992), 196-224.

[8] F. Bergeron and C. Reutenauer, Combinatorial resolution of systems of differential equations III: a special class of differentiably algebraic series. Europ. J. combinatorics, 11 (1990) 501-12.

[9] M. Bousquet-Mélou, New enumerative results on two-dimensional directed animals, Disc. Math 180 (1998) 73-106.

[10] M. Bousquet-Mélou and A.J. Guttmann, Enumeration of threedimensional convex polygons, Annals of combinatorics, 1, (1997), 27-53.

[11] M. Bousquet-Mélou and A.J. Guttmann, Three-dimensional selfavoiding convex polygons, Phys. Rev. E, 55, (1997), R6323-6.

[12] C. H. Chang, The spontaneous magnetisation of a two-dimensional rectangular Ising model, Phys. Rev. 88, 1422, (1952).

[13] C. A. W. Citteur and P. W. Kasteleyn, Critical behaviour of Ising systems on extremely anisotropic lattices, Physica 62, 17-40, (1972), ibid. 64, 415-26, (1973), ibid. 66, 94-110 (1973), ibid. 68, 491-510 (1973).

[14] A. R. Conway, Private communication.

[15] A. R. Conway, R. Brak and A. J. Guttmann, Directed animals on two-dimensional lattices, J. Phys. A. 26, (1993),3085-91.

[16] A. R. Conway and A. J. Guttmann, Square Lattice Self-Avoiding Walks and Corrections to Scaling, Phys. Rev. Letts. 77, 5284, (1996).

[17] A. R. Conway, A. J. Guttmann and M. P. Delest, The number of three-choice polygons, Mathl. Comput. Modelling 26 51-8, (1997)

[18] T. de Neef and I. G. Enting, Series expansions from the finite lattice method, J. Phys. A: Math. Gen. 10, 801-5, (1977).

[19] D. Dhar, M. H. Phani and M. Barma, Enumeration of directed site animals on two-dimensional lattices, J. Phys. A. 15, (1982) L279-84.

[20] I. G. Enting and A. J. Guttmann, In preparation.

[21] I. G. Enting, Series Expansions from the Finite Lattice Method, Nucl. Phys. B (Proc. Suppl.) 47, (1996) 180-7.

[22] M. L. Glasser, A. J. Guttmann, and W. P. Orrick, (In preparation).

[23] A. J. Guttmann and I. G. Enting, Solvability of some statistical mechanical systems, Phys. Rev. Letts. 76, (1996), 344-7.

[24] D. Hansel, J. M. Maillard, J. Oitmaa and M. J. Velgakis, Analytical Properties of the Anisotropic Cubic Ising Model, J. Stat. Phys, 48, (1987), 69-80.

[25] B. D. Hughes, Random walks and random environments: 1 Random walks, Clarendon Press (Oxford) 1995.

[26] M. T. Jaekel and J. M. Maillard, Symmetry relations in exactly soluble models, J. Phys. A: Math. Gen. 15, 1309-25, (1982).

[27] M. T. Jaekel and J. M. Maillard, Inverse functional relation on the Potts model, J. Phys. A: Math. Gen. 15, 2241-57, (1982).

[28] M. T. Jaekel and J. M. Maillard, A disorder solution for a cubic Ising model, J. Phys. A: Math. Gen. 18, 641-51, (1985).

[29] M. T. Jaekel and J. M. Maillard, A criterion for disorder solutions of spin models, J. Phys. A: Math. Gen. 18, 1229-38, (1985).

[30] I. Jensen and A. J. Guttmann, In preparation.

[31] L. Lipshitz D-finite power series. J. Algebra, 122 (1989) 353-73.

[32] N. Madras and G. Slade, The self-avoiding walk, Birkhäuser (Boston) (1993)

[33] J. M. Maillard, The inversion relation, J. Physique 46, 32a-419, (1984).

[34] L Onsager, Crystal statistics I: A two-dimensional model with an orderdisorder transition, Phys. Rev. 65 (1944) 117-49.

[35] A. Rechnitzer and W. P. Orrick, Private communication.

[36] A. L. Owczarek and T. Prellberg, Exact solution of the discrete (1+1)- dimensional SOS model with field and surface interactions, J. Stat. Phys. 70, 1175-94, (1993).

[37] A. Rechnitzer, W. P. Orrick and A. J. Guttmann, (In preparation).

[38] R. P. Stanley, Differentiably finite power series, Europ. J. Comb., 1, (1980), 175-88.

[39] Y. G. Stroganov, A new calculation method for partition functions in some lattice models, Phys. Letts. 74A, 116-8 (1979).

[40] I. Syozi and S. Naya, Symmetrical properties of two-dimensional Ising lattices, Prog. Theor. Phys. 24 829, (1960).

[41] C J Thompson, Classical Equilibrium Statstical Mechanics (1988) (The Clarendon Press, Oxford).

[42] H. Tsukahara and T. Inami, Test of Guttmann and Enting’s conjecture in the Eight Vertex Model, J. Phys. Soc. Japan 67, 1067-1070, (1998).

[43] A. Wiles, Modular Elliptic Curves and Fermat’s Last Theorem, Annals of Maths., 142, (1995), 443-51.

[44] T. T. Wu, B. M. McCoy, C. A. Tracy and E. Barouch, Spin-spin correlation functions for the two-dimensional Ising model: Exact theory in the scaling region. Phys. Rev. 13B, 316-74, (1976).