A Simple Rederivation of Onsager’s Solution of the 2D Ising Model Using Experimental Mathematics

I n this case study, we illustrate the great potential of experimental mathematics and symbolic computation by rederiving, ab initio, Onsager’s celebrated solution of the two-dimensional Ising model in zero magnetic field. Onsager’s derivation is extremely complicated, as are all the subsequent proofs. Unlike Onsager’s, our derivation is not rigorous, yet it is absolutely certain (and would be even if Onsager had not already derived it), and should have been acceptable to physicists who do not share mathematicians’ fanatical (and often misplaced) insistence on rigor.


MANUEL KAUERS AND DORON ZEILBERGER
I I n this case study, we illustrate the great potential of experimental mathematics and symbolic computation by rederiving, ab initio, Onsager's celebrated solution of the two-dimensional Ising model in zero magnetic field. Onsager's derivation is extremely complicated, as are all the subsequent proofs. Unlike Onsager's, our derivation is not rigorous, yet it is absolutely certain (and would be even if Onsager had not already derived it), and should have been acceptable to physicists who do not share mathematicians' fanatical (and often misplaced) insistence on rigor. (We make the convention that if i is n 1 , then i þ 1 ¼ 1, and if j ¼ n 2 , then j þ 1 ¼ 1.) DEFINITION 2. Let Mðn 1 ; n 2 Þ be the set of n 1 Â n 2 matrices whose entries are either 1 or À1 (of course, there are 2 n1n2 such matrices). The Laurent polynomial P n 1 ;n 2 ðx; yÞ is defined as follows:  DEFINITION 5. Let Mðn 1 ; n 2 ; n 3 Þ be the set of n 1 Â n 2 Â n 3 three-dimensional arrays whose entries are either 1 or À1 (of course, there are 2 n1n2n3 such arrays

Onsager's Solution
In 1944, Lars Onsager famously derived, and rigorously proved, the special case y ¼ 1 of Exercise 1.
Onsager's Explicit Formula for the Zero-Field 2D Ising Model. Let r : Onsager's proof [5], and all subsequent proofs, are very complicated. We will soon show how this formula could have been naturally derived, way back in 1941, if they had had the software and hardware that we have today (and even, probably, thirty years ago). Unlike Onsager's derivation, which is fully rigorous, ours is not. So from a strictly (currently mainstream) mathematical viewpoint, it would have been considered ''only'' a conjecture had it been done before Onsager's rigorous derivation. But this conjecture would have been so plausible that it would have been wholeheartedly accepted by the theoretical physics community.
Onsger's elegant solution involves an infinite series, which entails taking a limit. The definition of the function f(x) also involves taking a limit (namely of logðP n;n ðxÞÞ n 2 as n ! 1). Why is the former limit better than the latter?
Indeed, the notion of ''explicit,'' or ''closed form,'' is vague and cultural. In ancient Greece, a geometric construction was acceptable only if it used straightedge and compass. In algebra, for a long time, a solution was acceptable only if it could be expressed in terms of the four elementary operations and root extractions. In enumerative combinatorics, a solution was (and sometimes still is) considered to be in closed form only if it is a product and/or quotient of factorials. And there are many other examples.
In a famous position paper [7], Herb Wilf tackled this problem in combinatorics. He was inspired to write it when he was asked to referee a paper containing a ''formula'' for a certain quantity. It turned out that computing the quantity via the formula took much longer than using the definition. Inspired by the-at the timenew paradigm of ''computational complexity,'' he suggested that an ''answer'' is an efficient algorithm to compute the quantity in question.
How would we compute f(x), using the definition, for a specific ''numeric'' x? We can, in principle, compute the sequence of Laurent polynomials P n;n ðxÞ directly, for, say, n 30, get the finite sequence of numbers flog P n;n ðxÞ= n 2 g 30 n¼1 , see whether they get closer and closer, and estimate the limit. Alas, computing P n;n ðxÞ by brute force involves adding up 2 n 2 terms, each of which takes Oðn 2 Þ operations to compute. This is hopeless! Also, to be fully rigorous, one has to be able to find a priori bounds for the error, and for each find (rigorously) an n such that jf ðxÞ À logð P n;n ðxÞÞ=n 2 j\ for n ! n . This is truly hopeless.
On the other hand, using elementary calculus, Onsager's solution enables us to compute f(x) very fast, to any desired accuracy.
More importantly, physicists do not really care about the explicit form of f(x) (or more generally, the still wide open f(x, y), and g(x)); they want to know the exact locations of the singularities (critical points), which describe at what values of x (and hence at what temperatures) a phase transition occurs, e.g., at what temperature water boils or freezes. Even more importantly, they care about the nature of the singularities, in other words, how water boils rather than at what temperature (which depends, e.g., on pressure). From Onsager's solution, one can easily find, using elementary calculus, the location, and nature, of the singularity of G(z), and hence of f(x). It is impossible to extract this information directly from the definition.
This motivation may be interesting, but it is irrelevant to us. All we want is to answer Exercise 1 in the special case y ¼ 1, with as little effort as possible and making full use of the computer. We require only elementary calculus and very elementary matrix algebra. We don't even use eigenvalues!

Recommended Reading
Even though it is irrelevant to our story, for those readers who wish to know the context and background, we strongly recommend Barry Cipra's very lucid and very engaging introduction to the Ising model [1]. We also recommend the excellent books [6] and [9].

Symbol-Crunching
Of course, it would be nice to find an expression for f(x) in terms of the symbol x. Computing P n;n ðxÞ for any specific n is a finite (albeit huge) computation, involving summing 2 n 2 monomials, so we can't go very far. But let's assume that we live in an ideal world, or that quantum computing has become a reality. Then computing P n 1 ;n 2 ðxÞ, and in particular P n;n ðxÞ, being finite, is always possible. The first, very natural, step, already proposed in 1941, which was motivated by the combinatorial approach (see below and [6, Chapter 6, eq. (1.9)], where we replace x 2 by x), is to write where w ¼ xÀ1 xþ1 . It follows from a simple combinatorial argument that Z n1;n2 ðwÞ is a polynomial in w, of degree n 1 n 2 .
Taking logarithms and dividing by n 1 n 2 , we get log P n 1 ;n 2 ðxÞ n 1 n 2 ¼ À log 2 þ logðx À1 þ 2 þ xÞ þ log Z n 1 ;n 2 ðwÞ n 1 n 2 : Using the fact (do it!) that x À1 þ 2 þ x ¼ 4 1Àw 2 , we get that f ðxÞ ¼ log 2 À logð1 À w 2 Þ þ lim n!1 log Z n;n ðwÞ n 2 ; where w ¼ xÀ1 xþ1 . So from now, all we need is to find It turns out (and it follows from elementary considerations) that the sequence 1 n 2 log Z n;n ðwÞ converges in the sense of formal power series. More precisely, for every positive integer r, the coefficient of w r in F(w) (our object of desire) coincides with that of 1 n 2 log Z n;n ðwÞ as soon as n [ r. So a natural experimental mathematics approach would be to try to find as many Taylor coefficients of F(w) as our computer allows and look for a pattern that would enable us to conjecture a closed-form expression for the Taylor coefficients of F(w), thereby determining F(w) and hence f(x).
In an ideal world, with an indefinitely large computer, this very naive approach would have succeeded. Alas, as it turned out, we would have needed to compute P n;n ðxÞ for n ¼ 96, and since 2 96 2 is such a big number, this very naive brute-force approach is doomed to failure in our tiny universe.

Using Transfer Matrices
A much more efficient approach to computing the Laurent polynomials P n1;n2 ðxÞ (and hence the polynomials Z n1;n2 ðwÞ) was suggested in the seminal paper of Kramers and Wannier [3]. That was also Onsager's starting point. It is easy to see (see [6, p. 118]) that for each n 1 , there are easily computed 2 n 1 Â 2 n 1 matrices, let's call them A n 1 ðxÞ, such that P n1;n2 ðxÞ ¼ traceA n1 ðxÞ n2 : With today's computers, it is possible to compute these for n 1 12 and as large an n 2 as desired. But once again, one can (still) not go very far.
In 1941, B. L. van der Waerden suggested an ingenious (very elementary) combinatorial approach, described beautifully in Barry Cipra's article [1] (see also [6,Chapter 6] and [9] for nice accounts). He observed that the coefficients of w in the polynomial Z n 1 ;n 2 ðwÞ have a nice combinatorial interpretation. Putting N ¼ n 1 n 2 , it turned out (and is very easy to see; see [6]) that for every positive integer r, the coefficient of w r in Z n 1 ;n 2 ðwÞ, let's call it p r , is the number of ''lattice polygons'' with r edges that can lie in an n 1 Â n 2 ''toroidal rectangle,'' i.e., the set f0; . . .; n 1 g Â f0; . . .; n 2 g with 0 identified with n 1 and n 2 respectively. A lattice polygon is a collection of edges such that every participating vertex has an even number (0, 2, or 4) of neighbors. It follows in particular that p r is zero if r is odd.
It also follows from elementary combinatorial considerations that for n 1 ; n 2 [ r, the coefficient p r is a certain polynomial in N [6, p. 150, eq. (1.17)], and hence may be written p r ðN Þ, and we can write Now it also follows from elementary considerations, already known in 1941, that once you take the log, divide by N ¼ n 1 n 2 , and take the limit, only the coefficients of N in these ''Ising polynomials'' survive, and that F ðwÞ ¼ lim n!1 logðZ n;n ðwÞÞ n 2 ¼ X 1 r¼1 a ð1Þ r w r : It remains to compute as many Ising polynomials p r ðN Þ as our computers will allow us, extract the coefficients a ð1Þ r of N, and hope to detect a pattern that will enable us to conjecture the general coefficient of F(w), and hence know f(x).

How to Compute the Combinatorial Ising Polynomials
The first thing that comes to mind, and works well for small r, is to actually look for the kind of lattice polygons that can show up; but as r gets larger, this gets out of hand. Rather than do the intricate combinatorics, we use the fact that P n 1 ;n 2 ðxÞ ¼ traceA n 1 ðxÞ n 2 , from which we can compute Ó 2018 The Author(s), Volume 41, Number 1, 2019 Z n1;n2 ðwÞ for n 1 12 (say) and n 2 as large as desired. For each individual coefficient of w r (r even), we output it for sufficiently many specific n 1 and n 2 , and then using undetermined coefficients or interpolation, we fit them into a polynomial (whose degree we know beforehand). In fact, it is possible to get p 2r ðN Þ by looking at n 1 ¼ r À 2, n 2 [ r, by excluding obvious polygons that belong to the ðr À 2Þ Â n 2 toroidal rectangle but are impossible for a larger rectangle.

The Ising Polynomials
Using this very naive approach (using only matrix multiplication and then taking the trace), our beloved computers came up with the following first 10 Ising polynomials (we were able to find quite a few more, but as we will soon see, the first ten polynomials suffice): Hence F(w) begins thus: However, these ten terms (and even forty of them) do not suffice to suggest a pattern.

Duality Saves the Day
Way back in 1941, in the seminal paper that we have already mentioned, Kramers and Wannier discovered the duality relation (see [1] for a lucid explanation) the duality relation can be written as or in a more symmetric form, It follows that a more natural, and hopefully more userfriendly, function to consider is f ðxÞ :¼ f ðxÞ À log x þ x À1 À Á ; and we have that f ðxÞ is unchanged under the involution It is natural to change from the variable w to one that is invariant under the interchange x $ x Ã . There are many possibilities. Obviously, in order to ensure the invariance, we can set z ¼ Rðx; x Ã Þ for any symmetric rational function R. We only need to ensure that when w is expressed as a series in z, this series has positive order, so that we are allowed to substitute it into F(w). Since F(w) has only even exponents, we may also prefer that the series w ¼ wðzÞ have only odd exponents in z, so that the substitution does not introduce odd exponents into F(w).
If we try a generic template (ansatz) for a symmetric rational function with numerator and denominator of degree at most two, where a i;j and b i;j are undetermined coefficients, we get a system of polynomial equations that can be easily solved using Gröbner bases. This gets translated into an equation relating z and w by eliminating x, using the fact that x ¼ 1þw 1Àw . The (computer-generated) result is an equation of the form where the dots stand for certain linear combinations of the undetermined coefficients, which we suppress here because of their size. In order to ensure that the solution for w of this equation is a series in z with odd exponents only, it suffices to force the coefficients of all terms w i z j with i þ j even to equal zero. This gives a linear system whose solution brings the equation down to ðw À 1Þwðw þ 1Þða 0;0 þ a 0;1 þ a 0;2 Þ þ ð1 þ w 2 Þ 2 zðb 0;0 À b 1;0 þ b 2;0 Þ ¼ 0: z ¼ a 0;0 þ a 1;0 ðx þ x Ã Þ þ a 0;1 xx Ã þ a 2;0 ðx þ x Ã Þ 2 þ a 1;1 ðx þ x Ã Þxx Ã þ a 0;2 ðxx Ã Þ 2 b 0;0 þ b 1;0 ðx þ x Ã Þ þ b 0;1 xx Ã þ b 2;0 ðx þ x Ã Þ 2 þ b 1;1 ðx þ x Ã Þxx Ã þ b 0;2 ðxx Ã Þ 2 ; with unknown integer coefficients a i ; b i to be determined. So for a specific point x, the task is to find a so-called integer relation of the real numbers mðxÞ; . . .; x 10 mðxÞ; m 0 ðxÞ; . . .; x 10 m 0 ðxÞ. There are well-known algorithms for finding such relations [2,4]. In order to recover the relation from the values at a single point x, we would need to compute these values to a rather high precision, which is not an easy thing to do. We can get along with less precision using several evaluation points and searching for a simultaneous integer relation of the numbers mðxÞ; . . .; x 10 mðxÞ; m 0 ðxÞ; . . .; x 10 m 0 ðxÞ, for several x. It turns out that by using enough evaluation points, we need only about six decimal digits of accuracy of m(x) and m 0 ðxÞ for each of these points in order to establish a convincing guess. Unfortunately, this is a still bit more than what we were able to obtain by a direct computation via transfer matrices.

Supporting Software
For Maple and C programs, as well as output files, please visit the web page http://sites.math.rutgers.edu/*zeilberg/ mamarim/mamarimhtml/onsager.html.