A Simple Re-Derivation of Onsager’s Solution of the 2D Ising Model using Experimental Mathematics

Manuel KAUERS and Doron ZEILBERGER

Abstract: In this case study, we illustrate the great potential of experimental mathematics and symbolic computation, by rederiving, ab initio, Onsager’s celebrated solution of the twodimensional 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 (even if Onsager did not do it before), and should have been acceptable to physicists who do not share mathematicians’ fanatical (and often misplaced) insistence on rigor.

Two Warm-Up Exercises

Definition 1: For an n_{1} × n_{2} matrix, M = (m_{i,j}), and any positive real numbers x and y:

weight(M) (x, y) := x^{^{\frac{1}{2}} (∑_{i,j} m_{i,j} m_{i+1,j} + m_{i,j} m_{i,j+1})} · y^{∑_{i,j} m_{i,j}} .

(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^{n_{1} n_{2}} such matrices). The Laurent polynomial P_{n_{1},n_{2}}(x, y) is defined as follows.

P_{n_{1},n_{2}}(x, y) := ∑_{M∈M(n_{1},n_{2})} weight(M) (x, y) .

Definition 3: For x, y positive real numbers:

f(x, y) := lim_{n→∞} \frac{log P_{n,n}(x, y)}{n^{2}} .

Exercise 1: Find an explicit, closed-form, expression for f(x, y).

Definition 1a: For an n_{1} × n_{2} × n_{3} three-dimensional array, M = (m_{i,j,k}), and a positive real number x,

weight′(M) (x) := x^{^{\frac{1}{2}} ( ∑_{i,j,k} m_{i,j,k} m_{i+1,j,k} + m_{i,j,k} m_{i,j+1,k} + m_{i,j,k} m_{i,j,k+1} )} .

Definition 2a: 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^{n_{1} n_{2} n_{3}} such arrays), define the Laurent polynomial in x, by

Q_{n_{1},n_{2},n_{3}}(x) := ∑_{M∈M(n_{1},n_{2},n_{3})} weight′(M) (x) .

Definition 3a: For x, a positive real number,

g(x) := lim_{n→∞} \frac{log Q_{n,n,n}(x)}{n^{3}} .

Exercise 1a: Find an explicit, closed-form, expression for g(x).

We hope, dear readers, that you will spend some time trying to solve these two exercises, but please do not spend too much time! While we do know that the limits exist [vH], evaluating them explicitly has been an open problem for almost eighty years, and in spite of many attempts by the best minds in mathematical physics, both ‘excercises’ are still wide open.

Exercise 1 is called “solving the two-dimensional Ising model with magnetic field”, while Exercise 1a is called “solving the three-dimensional Ising model in zero magnetic field”. Let us quote Ken Wilson, who got the Physics Nobel prize in 1982 for seminal (non-rigorous!) work on questions related to these two ‘exercises’.

“When I entered graduate school I had carried out the instructions given to me by my father

[notable chemist E. Bright Wilson, who co-authored, with Linus Pauling, the classic Introduction to Quantum

Mechanics] and had knocked on both Murray Gell-Mann’s and Feynman’s doors and asked them what they were currently doing. Murray wrote down the partition function for the three-dimensional Ising model and said that it would be nice if I could solve it. Feynman’s answer was ‘nothing’ . ” [Quoted in Julia Yeomans’ wonderful book [Y], p. 35 .]

Onsager’s Solution

In 1944, Lars Onsager famously derived, and rigorously proved, the special case of Exercise 1, when y = 1.

Onsager’s Explicit Formula For the Zero-Field 2D Ising Model: Let

G(z) := −\frac{1}{4} ∑_{r=1}^{∞} \binom{2r}{r}^{2} \frac{z^{2r}}{r} ,

then

f(x, 1) = ln(x + x^{−1}) + G ( \frac{x − x^{−1}}{(x + x^{−1})^{2}} ) .

Onsager’s proof [O], and all the subsequent proofs, are very complicated. We will soon show how this formula could have been naturally derived, way back in 1941, if they had the software and hardware that we have today (and even, probably, thirty years ago).

Unlike Onsager’s derivation, that is fully rigorous, ours is not. So from a strictly (currently mainstream) mathematical viewpoint, it would have been considered ‘only’ a conjecture, were it done before Onsager’s rigorous derivation. But this conjecture would have been so plausible that it would have been whole-heartedly accepted by the theoretical physics community.

What is an “Explicit” Answer

From now on we will write f(x) instead of f(x, 1), and P_{n_{1},n_{2}}(x) instead of P_{n_{1},n_{2}}(x, 1).

Onsger’s elegant solution involves an infinite series, that entails taking a limit. The definition of the function f(x) also involves taking a limit (namely of \frac{log(P_{n,n}(x))}{n^{2}} as n → ∞). 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 geometrical construction was acceptable only if it used ruler 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 cominatorics, a solution was (and sometimes still is) considered closed form only if it is a product and/or quotient of factorials, and there are many other examples.

In a famous position paper [W], 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 time—new 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, and get the finite sequence of numbers {log P_{n,n}(x)/n^{2}}_{n=1}^{30}, and 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 take 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 |f(x) − log( P_{n,n}(x) )/n^{2}| < ϵ 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, physicist 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 location of the singularities, (critical points) that describes at what value of x (and hence at what temperature) 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 (that depends, e.g. on pressure). From Onsager’s solution, one can easily find, using Calculus I, 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 only require 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 do wish to know the context and background, we strongly recommend Barry Cipra’s [C] very lucid and very engaging introduction to the Ising model. We also recommend the excellent books [T] and [Y].

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 became 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, that was motivated by the combinatorial approach (see later, and [T], Ch.6, Eq. 1.9, where we replace x^{2} by x) is to write

P_{n_{1},n_{2}}(x) = \frac{(x + 2 + x^{−1})^{n_{1} n_{2}}}{2^{n_{1} n_{2}}} Z_{n_{1},n_{2}}(w) , where w = \frac{x − 1}{x + 1} .

It follows from a simple combinatorial argument that Z_{n_{1},n_{2}}(w) is a polynomial in w, of degree n_{1} n_{2}.

Taking logarithms, and dividing by n_{1} n_{2}, we get

\frac{log P_{n_{1},n_{2}}(x)}{n_{1} n_{2}} = − log 2 + log(x^{−1} + 2 + x) + \frac{log Z_{n_{1},n_{2}}(w)}{n_{1} n_{2}} .

Using the fact (do it!) that x^{−1} + 2 + x = \frac{4}{1−w^{2}} we get that

f(x) = log 2 − log(1 − w^{2}) + lim_{n→∞} \frac{log Z_{n,n}(w)}{n^{2}} where w = \frac{x − 1}{x + 1}.

So from now, all we need is to find

F (w) := lim_{n→∞} \frac{log Z_{n,n}(w)}{n^{2}} .

Now, it turns out (and it follows from elementary considerations) that the sequence \frac{log Z_{n,n}(w)}{n^{2}} converges in the sense of ‘formal power series’. More precisely, for any positive integer, r, the coefficient of w^{r} in F (w) (our object of desire) coincides with that of \frac{log Z_{n,n}(w)}{n^{2}} as soon as n > r. So a natural experimental mathematics approach would be to try and find as many Taylor coefficients of F (w) as our computer would allow 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_{n_{1},n_{2}}(x) (and hence the polynomials Z_{n_{1},n_{2}}(w)), was suggested in the seminal paper of Kramers and Wannier [KW]. That was also Onsager’s starting point. It is easy to see (see [T], p. 118) that for each n_{1}, there are easily computed 2^{n_{1}} by 2^{n_{1}} matrices, let’s call them A_{n_{1}}(x) such that

P_{n_{1},n_{2}}(x) = Trace A_{n_{1}}(x)^{n_{2}} .

With today’s computers, it is possible to compute these for n_{1} ≤ 12 and as large as 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 [C] (see also Chapters 6 of [T] and [Y] 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 [T]) that for any 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} by n_{2} ‘torodial rectangle’, i.e. the set {0, …, n_{1}} × {0, …, n_{2}} 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} > r, n_{2} > r, the coefficient p_{r} is a certain polynomial in N ([T], p. 150, Eq. (1.17)), and hence may be written p_{r}(N), and we can write:

p_{r}(N) = Na_{r}^{(1)} + N^{2}a_{r}^{(2)} + … + N^{m}a_{r}^{(m)} .

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→∞} \frac{log( Z_{n,n}(w) )}{n^{2}} = ∑_{r=1}^{∞} a_{r}^{(1)}w^{r} .

It remains to compute as many Ising polynomials, p_{r}(N), as our computers will allow us, extract the coefficients a_{r}^{(1)} of N, and hope to detect a pattern, to 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) = Trace A_{n_{1}}(x)^{n_{2}}, from which we can compute Z_{n_{1},n_{2}}(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} torodial rectangle but are impossible for a larger rectangle.

The Ising Polynomials

By using this very naive approach (only using 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).

p_{2}(N) = 0, p_{4}(N) = N, p_{6}(N) = 2 N, p_{8}(N) = \frac{1}{2} N (9 + N) , p_{10}(N) = 2 N (6 + N) ,

p_{12}(N) = \frac{1}{6} N (7 + N) (32 + N) , p_{14}(N) = N (130 + 21 N + N^{2}) ,

p_{16}(N) = \frac{1}{24} N (11766 + 1715 N + 102 N^{2} + N^{3}) , p_{18}(N) = \frac{1}{3} N (5876 + 776 N + 49 N^{2} + N^{3}) ,

p_{20}(N) = \frac{1}{120} N (980904 + 118830 N + 7415 N^{2} + 210 N^{3} + N^{4}) .

Extracting the coefficients of N, we get

0, 1, 2, \frac{9}{2}, 12, \frac{112}{3}, 130, \frac{1961}{4}, \frac{5876}{3}, \frac{40871}{5} .

Hence F (w) starts with

F (w) = w^{4} + 2 w^{6} + \frac{9}{2} w^{8} + 12 w^{10} + \frac{112}{3} w^{12} + 130 w^{14} + \frac{1961}{4} w^{16} + \frac{5876}{3} w^{18} + \frac{40871}{5} w^{20} + · · · .

However, these ten terms (and even forty of them) do not suffice to guess a pattern.

Duality Saves the Day

Way back in 1941, in the seminal paper of Kramers and Wannier, that we have already mentioned, they discovered the duality relation (see [C] for a lucid explanation)

f ( \frac{x + 1}{x − 1} ) = f(x) − log ( \frac{x − x^{−1}}{2} ) .

Letting

x^{∗} = \frac{x + 1}{x − 1} ,

the duality relation can be written as

f(x^{∗}) = f(x) − log ( \frac{x − x^{−1}}{2} ) ,

or in a more symmetric form

f(x) − log(x + x^{−1}) = f(x^{∗}) − log(x^{∗} + (x^{∗})^{−1}) .

It follows that a more natural, and hopefully user-friendly, function to consider is

f̄(x) := f(x) − log(x + x^{−1}) ,

and we have that f̄(x) is unchanged under the involution x ↔ x^{∗},

f̄(x^{∗}) = f̄(x) .

It is natural to change from the variable w to one that is invariant under the change 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) has 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,

z = \frac{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}}

where a_{i,j} and b_{i,j} are undetermined coefficients, we get a system of polynomial equations that can be easily solved using so-called Gröbner bases. This gets translated into an equation relating z and w by eliminating x, using the fact that x = \frac{1+w}{1−w}. The (computer-generated) result is an equation of the form

(…)+(…)w +(…)w^{2} +(…)w^{3} +(…)w^{4} +(…)z +(…)wz +(…)w^{2}z +(…)w^{3}z +(…)w^{4}z = 0 ,

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 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.

This suggests the choice

z = \frac{cw(1 − w^{2})}{(1 + w^{2})^{2}}, or w = \frac{z}{c} + \frac{3z^{3}}{c^{3}} + \frac{22z^{5}}{c^{5}} + \frac{211z^{7}}{c^{7}} + \frac{2306z^{9}}{c^{9}} + ⋯,

for some nonzero constant c. The value of c is not important. We take c = 2.

Let f̄(x), in terms of w, be written F̄(w), then (since x + x^{−1} = \frac{2(1+w^{2})}{(1−w^{2})}; note that x = \frac{1+w}{1−w})

F̄(w) := f(x) − log(x + x^{−1}) = − log(1 − w^{2}) + F (w) + log 2 − log ( \frac{2(1 + w^{2})}{1 − w^{2}} )

= − log(1 + w^{2}) + ∑_{r=0}^{∞} a_{r}^{(1)}w^{r} ,

giving

F̄(w) = −w^{2} + \frac{3}{2} w^{4} + \frac{5}{3} w^{6} + \frac{19}{4} w^{8} + \frac{59}{5} w^{10} + \frac{75}{2} w^{12}

+\frac{909}{7} w^{14} + \frac{3923}{8} w^{16} + \frac{17627}{9} w^{18} + \frac{81743}{10} w^{20} + O(w^{22}) .

Changing the variable to z, and renaming F̄(w) to G(z), we get

G(z) = −\frac{1}{4} z^{2} − \frac{9}{32} z^{4} − \frac{25}{48} z^{6} − \frac{1225}{1024} z^{8} − \frac{3969}{1280} z^{10} − \frac{17787}{2048} z^{12}

−\frac{184041}{7168} z^{14} − \frac{41409225}{524288} z^{16} − \frac{147744025}{589824} z^{18} − \frac{2133423721}{2621440} z^{20} + O(z^{22}) .

The first ten terms of the sequence of coefficients, let’s call them {b_{2r}}_{r=1}^{10}

−\frac{1}{4}, −\frac{9}{32}, −\frac{25}{48}, −\frac{1225}{1024}, −\frac{3969}{1280}, −\frac{17787}{2048}, −\frac{184041}{7168}, −\frac{41409225}{524288}, −\frac{147744025}{589824}, −\frac{2133423721}{2621440}, …

factorize nicely, which indicates that the series might be hypergeometric, i.e. the ratio of consecutive terms is a rational function of r. This is good news, since the famous sine and cosine functions and many other functions that come up in physics (e.g. the Hermite and Legendre polynomials that are so important in Quantum Mechanics) and elsewhere are hypergeometric.

By setting up an ‘ansatz’

\frac{b_{2r+2}}{b_{2r}} = \frac{a_{0} + a_{1}r + a_{2}r^{2} + a_{3}r^{3}}{b_{0} + b_{1}r + b_{2}r^{2} + r^{3}} ,

plugging-in the known values for 1 ≤ r ≤ 9, and simplifying, we get a system of nine linear equations with the seven unknowns a_{0}, a_{1}, a_{2}, a_{3}, b_{0}, b_{1}, b_{2}. If you take a random such system it is most likely unsolvable. If the computer finds a solution, it is great news. What is true for the first 9 values is probably true for ever.

This means that the sequence of ratios b_{2r+2}/b_{2r} probably matches a rational function in r. Given the ratios {b_{2r+2}/b_{2r}}_{r=1}^{9}, the computer immediately established that

\frac{b_{2r+2}}{b_{2r}} = \frac{r(2r + 1)^{2}}{(r + 1)^{3}} ,

for 1 ≤ r ≤ 9, and if true for all r, this would imply the closed-form expression, for the coefficients

b_{2r} = −\frac{\binom{2r}{r}2}{r4^{r+1}} .

Since we can (nowadays!) easily extend the sequence b_{2r} up to (at least) sixteen terms, and this ‘guess’ indeed continued to hold, this makes it virtually certain that the guess is correct. Combining everything, we derived, ab initio, by pure guessing (and very elementary and natural reasoning), Onsager’s formidable formula.

What’s next?

Now that we have rediscovered Onsager’s explicit formula for f(x) = f(x, 1), a natural next step towards the general case f(x, y) is to determine an explicit expression for m(x) = \frac{d}{dy} f(x, y)|_{y=1}, i.e., the next term in the Taylor series expansion of f(x, y) with respect to y at y = 1. Physicists call this the “spontaneous magnetization”.

Using transfer matrices, as before, it is easy to compute the first few terms of m(x) as a series in x (or w, or z), and we don’t even need a computer to guess an explicit expression for them: they all are zero. But that’s just a part of the story.

Onsager observed that m(x) is only zero for x < 1 + \sqrt{2}, while for x ≥ 1 + \sqrt{2}, it is equal to

(\frac{(x^{2} + 1)^{2}(x^{2} − 2x − 1)(x^{2} + 2x − 1)}{(x − 1)^{4}(x + 1)^{4}})^{1/8} .

According to Thompson ([T], p. 135), this expression “was first derived by Onsager in the middle of the 1940s, but in true Onsager fashion he has not to this day published his derivation”.

We don’t know how he found this expression, but here is one way one could search for it, using experimental mathematics. For specific numbers x, y, we can compute numerical approximations of f(x, y) using the original definition (Def. 3 above). For example, taking f(x, y) ≈ log P_{n,n}(x, y)/n^{2} with n ≈ 20 gives several correct digits at a reasonable computational cost. From the numerical estimates of f(x, y) for various points x, y, we can obtain numerical estimates for m(x) and m′(x), for various points x.

The idea is to fit a differential equation against this numeric data. Suppose we suspect a differential equation of the form

(a_{0} + a_{1}x + ⋯ + a_{10}x^{10}) m(x) + (b_{0} + b_{1}x + ⋯ + b_{10}x^{10}) m′(x) = 0,

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′(x), …, x^{10}m′(x). There are well-known algorithms for finding such relations [FB, LLL].

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 by using several evaluation points and searching for a simultaneous integer relation of the numbers m(x), …, x^{10}m(x), m′(x), …, x^{10}m′(x), for several x. It turns out that by using enough evaluation points, we just need about 6 decimal digits of accuracy of m(x) and m′(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 .

Acknowledgement. We thank the anonymous referee for many insightful remarks, that improved our paper considerably. References

[C] Barry A. Cipra, An Introduction to the Ising Model, The American Mathematical Monthly 94 (1987), 937–959.

Available from http://www.yaroslavvb.com/papers/cipra-introduction.pdf (accessed May

9, 2018).

[FB] H.R.P. Ferguson and D.H. Bailey, A Polynomial Time, Numerically Stable Integer Relation Algorithm. RNR Techn. Rept. RNR-91-032, Jul. 14, 1992.

[KW] H.A. Kramers and G.H. Wannier, Statistics of the Two-Dimensional Ferromagnet. Part I, Physical Review 60 (1941), 252-262.

Available from http://sites.math.rutgers.edu/~zeilberg/akherim/KW1941.pdf (accessed May

9, 2018).

[LLL] A.K. Lenstra, H.W. Lenstra, and L. Lovasz, Factoring Polynomials with Rational Coefficients. Math. Ann. 261, 515-534, 1982.

[O] Lars Onsager, Crystal Statistics, I. A two-dimensional model with an order-disorder transition, Physical Reviews 65(1944), 117–149.

Available from http://www.phys.ens.fr/~langlois/Onsager43.pdf (accessed May 9, 2018).

[T] Colin J. Thompson, “Mathematical Statistical Physics”, MacMillan, 1972.

[W] Herbert S. Wilf, What is an Answer?, The American Mathematical Monthly 89 (1982), 289– 292.

[vH] L. van Hove, Sur L’intgrale de Configuration Pour Les Systmes De Particules Une Dimension, Physica 16 (1950), 137–143.

[Y] J.M. Yeomans, “Statistical Mechanics of Phase Transition”, Oxford University Press, 1992.

Manuel Kauers, Institute of Algebra, J. Kepler University Linz, Altenbergerstrasse 69, 4040 Linz, Austria.

manuel dot kauers at jku dot at .

Doron Zeilberger, Department of Mathematics, Rutgers University (New Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA.

Email: DoronZeil at gmail dot com .

First written: May 22, 2018; This version: Sept. 04, 2018.