%begin macros
\parskip=8pt
\magnification=\magstephalf
\parindent=0pt
\overfullrule=0in
%end macros
\rm
 
\bf 
A Computational Approach To Finding Exact Solutions Of Ising-Type
Problems\rm 
 
{\it Doron Zeilberger}
\footnote{ *}{
Department of Mathematics and Computer Science, Drexel University,
Philadelphia, PA 19104. Supported in part by the NSF. Written 1989.
}
 
{\bf Never published before. Exclusive for DZ's website}
 
{\bf 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.
 
{\it PACS numbers:  05.50+q   64.60.Cn   02.60.+y   02.70+d.}
 
Computational methods are traditionally associated with the attempt
to find {\it 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 {\it computational} method for finding 
{\it exact} solutions is proposed. 
It is based on the realization that
known solutions of two-dimensional models belong to the 
class of so-called {\it D-finite} functions. 
A function $f(x)$ is {\it D-finite}
if it is a solution of a linear differential equation with \it 
polynomial
coefficients\rm 
 
$$
 a_0 (x) f(x) + a_1 (x) f'(x) + ... + a_n (x) f^{(n)} (x) = b(x) \, \, \,, 
\eqno(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 {\it 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 $\rho$,
that is obtained by substituting $f(x)= (x- x_0 )^{\rho}$
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
{\it Pade approximation}, 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 D-finite. Then
one can still get "D-finite approximations" as follows.
Compute the $( N_1 , ... , N_n ; M)$ \it 
D-finite
approximant\rm 
 as 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 {\it rational} function $R(x)$.
 
A function $f( x_1 , ... , x_n )$ is {\it holonomic} (or
{\it 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 
 
$$
 \int_{a}^{b} f( x_1 , ... , x_n ) d x_n  
\eqno(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
 
$$
 \Lambda ( v_1 , ... , v_r ) = \int_{0}^{2 \pi } 
\int_{0}^{ 2 \pi } lnP( e^{v_1} , ... , e^{v_r} , e^{ i \theta_1 }, e^{ i \theta_2 } )  d \theta_1 d \theta_2 
\eqno(3)
$$       
 
where $P$ is a polynomial. Now making the change of variables
$ z_1 = e^{ i \theta_1 }$,$ z_2 = e^{ i \theta_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 $\Lambda$ 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 $ \psi /kT$, since $U= J(1- w^2 ) {{{d}\over{dw}}} (- \psi /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} - 4 w^{10} +35 w^8 -56 w^6 + 15 w^4 - 20 w^2 -3)V(w) - 
$$
 
$$
  (1+ w^2 ) ( w^2 +2w-1)(-1-2w+ w^2 )(3 w^8 - 10 w^6 + 20 w^4 + 2 w^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) \equiv 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 {\it P-recursive} sequences. A sequence
$a(m)$ is {\it 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) \, \, \,, 
\eqno(4)
$$         
 
where $p_0 , ... , p_L , and r$ are certain polynomials
in $m$. One of the simplest "P-recursive" sequences are the binomial
coefficients $ {{n} \choose {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 {\it q-binomial coefficients}
$ {{n} \choose {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 {\it 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 two-dimensional
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 P-recursive, 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 \beta H)$.
 
What about higher dimensions?
It is well known that the solution of the one-dimensional Ising
model is an {\it 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 
$\infty times \infty \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 \rightarrow 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.
 
{\bf References}
 
[A] G. Andrews, "The Theory of Partitions", Addison-Wesley, Reading,
1976.
 
[Ba] G. A. Baker, Jr., {\it Phys. Rev.} {\bf 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) \equiv 0$), in (1).
 
[Bax] R. J. Baxter, {\it J. Phys. A} {\bf 13} (1980), L61.
 
[Be]  I. N. Bernstein, {\it Functional Analysis and Its applications}, 
{\bf 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, {\it Adv. Phys.} {\bf 9}(1960),150.
 
[G1] I. M. Gessel, {\it Symmetric functions and P-recursiveness},
J. Comb. Th., Ser. A, in press.
 
[G2] I. M. Gessel, {\it Enumerative applications of symmetric functions},
Proc. of the $17^{th}$ Lotharingien Seminar, IRMA, University of
Strasbourg, 1988.
 
[GJ] I. P. Goulden, D. M. Jackson, {\it Siam J. Alg. Disc. Meth.} {\bf 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, {\it Physica} {\bf 27}(1961), 1209.
 
[O] L. Onsager, Phys. Rev. {\bf 65}(1944),117.
 
[S] R. P. Stanley, Europ. J. Comb. {\bf 1}(1980), 175.
 
[TF] H. N. V. Temperley and M. E. Fisher, {\it Phil. Mag.} {\bf 6} (1961),
1061.
 
[T] C. J. Thompson, "Mathematical Statistical Mechanics", Princeton
University Press, 1979.
 
[Z] D. Zeilberger, \it 
A holonomic systems approach to special
function identities\rm 
, preprint.
 
\end
 
 
