A Computational Approach To Finding Exact Solutions Of Ising-Type Problems

Doron Zeilberger *

Never published before. Exclusive for DZ’s website

Abstract: A computational method for finding exact solutions of Ising type problems is proposed. It is shown that known solutions of two-dimensional Ising and dimer problems satisfy linear differential equations with polynomial coefficients, that could have been found computationally from a large enough series expansion. The knowledge of the differential equation is a legitimate form of expressing the exact solution since it enables one to find the critical point and exponent. This raises the hope that by using sufficiently large series expansions it would be possible to find, by computation, exact solutions of still open problems.

PACS numbers: 05.50+q 64.60.Cn 02.60.+y 02.70+d.

Computational methods are traditionally associated with the attempt to find approximate solutions to physical problems. On the other hand the search for ”exact” or ”analytic” solutions is considered an exclusively human, ”paper and pencil”, pursuit.

A computational method for finding exact solutions is proposed. It is based on the realization that known solutions of two-dimensional models belong to the class of so-called D-finite functions. A function f(x) is D-finite if it is a solution of a linear differential equation with polynomial coefficients

a_{0}(x)f(x) + a_{1}(x)f′(x) + ... + a_{n}(x)f^{(n)}(x) = b(x) ,  (1)

where a_{0}(x) , a_{1}(x) , ..., a_{n}(x), b(x) are polynomials.

Given a function that we know or conjecture is D-finite and given a certain number of terms, say L, of its Taylor expansion around x = 0, we can try and find the differential equation (1) it satisfies as follows. Guess an upper bound n for the order and an upper bound B for the degrees of the coefficients. If (n + 2)(B + 1) < L we can substitute the known portion of f(x) into (1), extract the coefficients of x^{i , i=0}...L − 1 and get L linear equations for the (n + 2)(B + 1) unknown coefficients of the a_{i}(x)′s and b(x). Since the number of equations exceeds the number of unknowns, there is no a priori reason to expect a non-zero solution. If we do find such a solution it would mean that beyond all reasonable doubt the sought-for function is indeed D-finite and we have found a linear differential equation satisfied by f(x). If we fail, we raise B or n. If (n + 2)(B + 1) gets to be bigger than L we have to increase L, i.e. obtain a larger series expansion.

Knowing a linear differential equation satisfied by the sought-for function is as good as having a ”closed form” solution. Once the differential equation is known one can easily determine exactly the critical point and the critical exponent[Bu]. The critical point is the smallest positive real root x_{0} of the equation a_{n}(x) = 0 (where a_{n}(x) is the leading coefficient of (1)), and the critical exponent is the solution of the ”indicial equation” in ρ, that is obtained by substituting f(x) = (x − x_{0})^{ρ} into (1) and equating to zero the lowest degree term. (In case of a double root we get a logarithmic singularity. )

The present method can be viewed as a generalization of the method of Pade approximation,

* Department of Mathematics and Computer Science, Drexel University, Philadelphia, PA 19104. Supported in part by the NSF. Written 1989. introduced to physics by Baker[Ba].

This naturally leads to the concept of ”D-finite approximation”. Suppose that one does not have enough terms to find the true differential equation, or that the function is not believed to be Dfinite. Then one can still get ”D-finite approximations” as follows. Compute the (N_{1}, ..., N_{n}; M) D-finite approximantas the differential equation (1), with degree(a_{i}) = N_{i}, such that if the first L := N_{1} + ...N_{n} + n + M + 1 terms in the series expansion of f are substituted into (1) one gets that all the coefficients up to x^{L−1} agree. It would be interesting to develop a general theory of these ”D-finite approximations” and to make numerical studies to see if, for example, the (3,3,3;0) D-finite approximant gives a better approximation for the critical point and exponent then the (6,6) Pade approximant.

Although the present method is non-rigorous, it might be possible to make it rigorous. By using the combinatorial interpretation or otherwise it is conceivable that it would be possible to prove ”from the outset”, using general theorems, that a certain function of interest is D-finite. This was done for large families of functions ([G1][G2][GJ][Z])

If in addition one can give a priori bounds for the order, and degrees of the coefficients, of the differential equation (1) that it satisfies, it would follow then that if the series expansion agrees to sufficiently many terms with the solution of the ”approximate” differential equation, we will have a rigorous proof that the ”approximate” differential equation is indeed exact.

The class of D-finite functions was introduced as such by Stanley[S]. It is the one-dimensional case of Bernstein’s class([Be][Bj])

It is well known and easy to see[Z][S] that the class of D-finite functions is an algebra: the sum and product of two D-finite functions is again D-finite. Furthermore it is also easy to see that if f(x) is D-finite so is f(R(x)), for any rational function R(x).

A function f(x_{1}, ..., x_{n}) is holonomic (or multi D-finite or ”belongs to the Bernstein class” [Z][Be][Bj] if, (possibly after a certain linear change of coordinates) it satisfies a system of ”ordinary” equations in each of its variables, such that the ”characteristic set” (i.e. the set of common zeros of the leading coefficients) is finite. It was proved[Z] (using ideas of Bernstein[Be] ) that the class of holonomic functions is an algebra, and that if f(x_{1}, ..., x_{n}) is holonomic in x_{1}, ..., x_{n}, then

∫_{a}^{b} f(x_{1}, ..., x_{n})dx_{n}  (2)

is holonomic in the surviving variables x_{1}, ..., x_{n−1}. By iterating (2), we get that the definite integral of a holonomic function f with respect to some of its variables is holonomic in its surviving variables.

The known solutions of the two-dimensional Ising problems for the square, triangular, and hexagonal lattices [O][D] and dimer problems[K][TF] (but not Baxter’s[Bax] hard hexagon) have the form

Λ(v_{1}, ..., v_{r}) = ∫_{0}^{2π} ∫_{0}^{2π} lnP (e^{v_{1}}, ..., e^{v_{r}}, e^{iθ_{1}}, e^{iθ_{2}})dθ_{1}dθ_{2}  (3)

where P is a polynomial. Now making the change of variables z_{1} = e^{iθ_{1}},z_{2} = e^{iθ_{2}}, and noting that lnP is holonomic for any polynomial P (because the derivative w.r.t to any of its variables is a rational function), we get that Λ above is holonomic, as a function of e^{v_{1}}, ..., e^{v_{r}}, and thus also as a function of w_{i} = tanh v_{i} , i = 1, ..., r. By setting nu_{1} = ... = nu_{r} = nu in (3) we see that the customary one-variable solutions of Ising-type problems are D-finite.

Baxter’s solution(Ba of the hard hexagon model has a completely different form. The solution kappa(z) is given parametrically, where both kappa and z are certain quotients of theta-like infinite products in an auxiliary variable q. Since the complete elliptic integral K(k), that is D-finite in k, can be so described[IK].

It is conceivable that Baxter’s solution is indeed D-finite. It should not be too hard to check this conjecture computationally, by using Baxter’s already known solution to extract the first 200 or so terms.

The claim that Onsager’s solution is D-finite will be now made concrete by exhibiting the differential equation it satisfies. For convenience we give the differential equation satisfied by V := wU − (1 + w^{2}), where (with the notation of ref. 14), U is the internal energy and w = tanh(v) = tanh(J/kT). ¿From the differential equation for V it is easy to get the differential equation for U and hence for the free energy per spin ψ/kT, since U = J(1 − w^{2})\frac{d}{dw} (−ψ/kT).

We used the expression for U given in formula (5.15), p.134 of Thompson’s book(Th in terms of the complete elliptic integral of the first kind K(k). Starting with the well known[BB] differential equation satisfied by K(k), and performing the appropriate change of independent and dependent variables, it was found, using MAPLE, that V (w) satisfies the following second-order linear differential equation:

4w(w^{12} − 4w^{10} + 35w^{8} − 56w^{6} + 15w^{4} − 20w^{2} − 3)V (w)−

(1 + w^{2})(w^{2} + 2w − 1)(−1 − 2w + w^{2})(3w^{8} − 10w^{6} + 20w^{4} + 2w^{2} + 1)V′(W)+

w(w − 1)(w + 1)(1 + w^{2})^{2}(w^{2} + 2w − 1)^{2}(−1 − 2w + w^{2})^{2}V′′(w) ≡ 0

Of course, the way we derived this equation is ”backwards”, by using the already known Onsager solution. Our point is, however, that Onsager’s solution does indeed satisfy a certain linear differential equation that is ”not too big”, and that could have easily been found computationally if one had the series expansion of the solution up to the w^{60} term.

What about non-zero external field? It is well known(St and easy to see that f(x) is D-finite if and only if its sequence of coefficients a(m) belongs to the class of so-called P-recursive sequences. A sequence a(m) is P-recursive if it satisfies a linear recurrence equation with polynomial coefficients:

p_{0}(m)a(m) + p_{1}(m)a(m + 1) + ... + p_{L}(m)a(m + L) = r(m) ,  (4)

where p_{0}, ..., p_{L}, andr are certain polynomials in m. One of the simplest ”P-recursive” sequences are the binomial coefficients \binom{n}{k}, which are P-recursive in both n and k, and that have an obvious interpretation as the number of paths with positive unit steps from (0, 0) to (k, n − k). It is well known[A] that if we attach a ”weight” of q(swept area) to every such path, then the generating function is given by the q-binomial coefficients \binom{n}{k}_{q} := (1 − q^{n})... (1 − q^{n−k+1})/(1 − q)...(1 − q^{k}). Now the q-binomial coefficients are the simplest example of a q-P-recursive sequence (in both k and n). A sequence a(m) of polynomials (or power series) in q is said to be q-P-recursive if it satisfies a linear recurrence equation with coefficients that are polynomials in q and q^{m}:

p_{0}(q, q^{m})a(m) + p_{1}(q, q^{m})a(m + 1) + ... + p_{L}(q, q^{m})a(m + L) = r(q, q^{m}) ,

It is known[T] that introducing area in van der Waerden’s combinatorial approach to the twodimensional Ising model would give the solution of the non-zero external field problem. Since we saw above that Onsager’s solution is D-finite, and thus its sequence of Mclaurin coefficients is Precursive, it is reasonable to conjecture that the still unknown solution of the more general problem will have a q-P-recursive sequence of coefficients and therefore the solution itself will be a so-called ”q-D-finite function”, i.e. satisfy a functional equation of the form f(x) + a_{1}(x, q)f(qx) + ... + a_{M}(x, q)f(q^{M}x, q) = 0. Hopefully the existing series expansion will lead to finding this equation computationally. Once we have a functional equation of the above form it is easy to find exactly the critical point and the critical exponent, that now depend on q = exp(−sβH).

What about higher dimensions? It is well known that the solution of the one-dimensional Ising model is an algebraic power series, i.e. there are polynomials P_{0}, ..., P_{L} such that P_{0}+...+P_{L}f^{L} = 0. The two-dimensional Ising model can be viewed as a limit of solutions of one dimensional problems on infinite strips of increasing height, but the degree of the polynomial equations that the solutions satisfy grows exponentially as the height of the strip. In the limit it ceases to be algebraic and turns into a D-finite function. By analogy, the solution of the ∞times∞ × L problem for every fixed L, is probably holonomic, but the corresponding differential equations may get to be of higher and higher order, and the limit, as L → infinity, which would be the solution of the three-dimensional problem, may belong to a yet to be discovered larger class of functions.

On the other hand, it is conceivable that the class of D-finite functions is already big enough that the solutions of the three- and higher-dimensional problems are D-finite. Gessel[G1] has recently shown that the number of graphs whose degrees belongs to a given finite set (for example 0,2,4,6) is P-recursive in both the number of vertices and the number of edges. Gessel’s results apply to general graphs, i.e. subgraphs of the complete graph. If it were true for subgraphs of the three-dimensional lattice, then it would mean that the solution of the three-dimensional Ising model is indeed D-finite and thus amenable to the present approach. Be that as it may, it would be interesting to find ”D-finite approximations” and compare their performance with Pade-approximations.

References

[A] G. Andrews, ”The Theory of Partitions”, Addison-Wesley, Reading, 1976.

[Ba] G. A. Baker, Jr., Phys. Rev. 124(1961), 768. Indeed finding the (N, N) Pade approximant of f′(x)/f(x) amounts to taking n = 1, B = N,L = 2(N + 1) (with b(x) ≡ 0), in (1).

[Bax] R. J. Baxter, J. Phys. A 13 (1980), L61.

[Be] I. N. Bernstein, Functional Analysis and Its applications, 5 (1971), 89.

[Bj] J. E. Bjork, ”Rings of Differential Operators”, North Holland, 1979.

[BB] J. M. Borwein and P. M. Borwein, ”Pi and the AGM”, J. Wiley, New York, 1987. p.9, formula (1.3.8).

[Bu] J. C. Burkill, ”The Theory of Ordinary Differential Equations”, Oliver and Boyd, Edinburgh and London (1962), Ch. 4.

[D] C. Domb, Adv. Phys. 9(1960),150.

[G1] I. M. Gessel, Symmetric functions and P-recursiveness, J. Comb. Th., Ser. A, in press.

[G2] I. M. Gessel, Enumerative applications of symmetric functions, Proc. of the 17^{th} Lotharingien Seminar, IRMA, University of Strasbourg, 1988.

[GJ] I. P. Goulden, D. M. Jackson, Siam J. Alg. Disc. Meth. 7 (1986), 179.

[IK] S. Iyanaga and Y. Kawada (editors), ”Encyclopedic Dictionary of Mathematics”, MIT Press, Cambridge, 1^{st} edition, 1980. Vol 2, App. A, Table 16.III, (1), (p. 1455).

[K] P. W. Kasteleyn, Physica 27(1961), 1209.

[O] L. Onsager, Phys. Rev. 65(1944),117.

[S] R. P. Stanley, Europ. J. Comb. 1(1980), 175.

[TF] H. N. V. Temperley and M. E. Fisher, Phil. Mag. 6 (1961), 1061.

[T] C. J. Thompson, ”Mathematical Statistical Mechanics”, Princeton University Press, 1979.

[Z] D. Zeilberger, A holonomic systems approach to special function identities, preprint.