SYMBOLIC COMPUTATION TO THE AID OF STATISTICAL MECHANICS

MANUEL KAUERS^{1} AND DORON ZEILBERGER

Abstract. We suggest that symbolic computation is a powerful tool in theoretical physics, that is currently not as widely used as it should be. We illustrate it by a bit of “alternative history”, how the celebrated and extremely complicated Onsager solution could (at least in principle, but probably also in pratice) could have been conjectured, using symbolic computation. Hopefully, the same methodology would be helpful to make progress on till open problems, like the 2D Ising model with a magnetic field, and the 3D Ising model.

1. Prologue

Until about thirty years ago, physics consisted of two parts: theoretical and experimental. Then a new kid came to the block, computational physics, that used a computer to ‘pretend’ that it is nature, and perform many experiments that would be too expensive, and usually impossible to perform in the real world. Alas, computational physics, with a few pioneering exceptions, used numerical computations (mostly simulations). We believe that there is room for yet another ‘kid’, symbolic computations. Very impressive applications of Symbolic Computation to high energy physics are described, for example, in [1]. Here we hope to illustrate the usefullness of symbolic computation to statistical mechanics, by focusing, as a case study, on the celebrated Ising model (e.g. [5]).

2. A crash course in statistical mechanics

The general scenario of statistical mechanics consists of very many possible configurations, each with a certain energy. Not all configurations are created equal, and some are much more likely than others, as first realized by the Austrian scientific giant, Ludwig Boltzmann. The probability of a configuration C, whose energy is E(C) is proportional to

e^{\frac{−E(C)}{kT}} ,

where k is Boltzmann’s constant, and T is the absolute temperature. It follows that the higher the energy, the less likely it is to show up, the decay being exponential. In order to get macroscopic (thermodynamic) (average) information, we form the famous Gibbs Partition function

Z = ∑_{C∈configurations} e^{\frac{−E(C)}{kT}} ,

and the ‘average’ energy (that physicists call internal energy) is obtained by taking the logarithmic derivative. Other quantities of interest are obtained by taking higher order derivatives. (Note that in taking the logarithmic derivative, Z shows up in the denominator, and this makes sense, in order to normalize the sum of the probabilities to be 1). Statistical mechanicians are fond of toy models that ‘simplify’ nature, and try to preserve the qualitative features of the real world. Unfortunately, even these toy models are extremely difficult, and with rare exceptions, intractable (or at least wide open), and people resort to approximations, using numerical methods, sophisticated brilliant (but non-rigorous) renormalization theory calculations, and most often, simulations (aka “Monte Carlo methods”).

Date: Feb. 15, 2018.

1 Partially supported by the Austrian FWF grant Y464-N18.

In this article we will show how symbolic computation may give a new approach to tackling these difficult models.

3. The grandaddy of all models: The (two dimensional) Ising Model

Fix positive integers m and n, and let M(m, n) be the set (of cardinality 2^{mn}) of all m×n matrices whose entries are −1 or 1 (called ‘spin down’ and ‘spin up’, respectively). Let

x = e^{\frac{J}{kT}} , y = e^{\frac{H}{kT}}

where J is the so-called coupling constant (and hence a (physical) constant) and H is another input variable (in addition to T) called the (strength of the) external magnetic field. The energy associated to a specific matrix M ∈ M(m, n) is defined as follows

E(M) := −J (∑_{i=1}^{m} ∑_{j=1}^{n} (M_{i,j}M_{i,j+1} + M_{i,j}M_{i,j+1})) − H (∑_{i=1}^{m} ∑_{j=1}^{n} M_{i,j}) .

Here we use the ‘periodic’ convention that M_{m+1,j} = M_{1,j} and M_{i,n+1} = M_{i,1}. It follows that the weight of a matrix is

weight(M) := exp(−\frac{E(M)}{kT}) = x^{∑{m \atop i=1} ∑{n \atop j=1}(M_{i,j}M_{i,j+1}+M_{i,j}M_{i,j+1})} · y^{∑_{i=1}^{m} ∑_{j=1}^{n} M_{i,j}} .

Our first (modest!) goal is just to compute, for as large m and n that our computers and human ingenuity would allow, the (Gibbs) Partition function (i.e. weight enumerator of the set M(m, n)).

Z_{m,n}(x, y) := ∑_{M∈M(m,n)} weight(M) ,

a certain Laurent polynomial in the variables x (of degree 2mn and low-degree ≥ −2mn) and y (of degree mn and low-degree −mn). Later on we would be especially interested in the ‘square case’ (m = n), and the semi-infinite strip case (n = ∞).

4. Transfer matrices

A direct [4] way to compute Z_{m,n}(x, y) is to simply construct the set M(m, n) consisting of all 2^{mn} whose entries are drawn from the set {−1, 1} matrices, compute the weight for each and every one of them, and add them all up. But there is a much more efficient way [5, 6], called the transfer-matrix method. Let’s associate for each integer i between 0 and 2^{m} − 1 the m-dimensional {−1, 1} column vector v(i) := (2b_{0} − 1, …, 2b_{m−1} − 1), where i = ∑_{j=0}^{m−1} b_{j}2^{j} is the binary representation of i. Define the 2^{m} × 2^{m} transfer matrix

T(m) := (T_{ij})_{0≤i,j≤2^{m}−1} ,

where the typical entry T_{ij} is as follows.

T_{ij} := x^{∑_{l=1}^{m} v(i)_{l}v(j)_{l}+∑_{l=1}^{m−1} v(j)_{l}v(j)_{l+1}+v(j)_{m}v(j)_{1}} · y^{∑_{l=1}^{m} v(j)_{l}} .

Using the transfer matrix, the desired polynomial Z_{m,n}(x, y) can be expressed very succinctly as

Z_{m,n}(x, y) = Tr [ T(m)^{n}] .

It is well-known (and easy to see!) that exponentiation of a matrix is really fast (poly-log time, using A^{2k} = (A^{k})^{2}, A^{2k+1} = A^{2k}A rather than A^{k} = A^{k−1}A), but for a numeric matrix. Here we have a symbolic matrix, and a naive approach, (once m gets larger) would be unfeasible. But with the method of homomorphic images, we can go pretty far.

5. The Functions that physicists care about

For the finite (torodial) m × n lattice, one is interested in the Free energy, F

F = −kT log Z ,

Internal Energy, U

U = kT^{2} \frac{∂}{∂T} log Z ,

and Specific heat

C = \frac{∂U}{∂T} .

Physicists are also very interested in the Magnetization, M

M = −\frac{∂F}{∂H} ,

and the Isothermal susceptibilty

χ_{T} = \frac{∂M}{∂H} .

At the “end of the day”, people take the limit as both m and n go to infinity, but one can already get a rough idea on what is going on for small values of m and n, as we hope to show soon.

6. The Data for Z_{m,n}(x, y)

For Z_{n,n}(x, y) for 1 ≤ n ≤ 9 see:

           http://www.math.rutgers.edu/~zeilberg/tokhniot/ising/Polys

Very impressively, Per Hakan Lundow[3], has already computed these for n ≤ 16 For Z_{m,n}(x, y) for 1 ≤ m ≤ 5 , 1 ≤ n ≤ 11 see:

           http://www.math.rutgers.edu/~zeilberg/tokhniot/ising/Polys1

7. The Semi-Infinite case

Physicists are not really interested in the polynomials Z_{m,n}(x, y) for their own sake but in the logarithm

log Z_{m,n} .

But, physicists (and chemists) that believe in the (fictional!) infinity, are really interested in the free energy F(x, y)

F := −kT lim_{m→∞} lim_{n→∞}(mn)^{−1} log Z_{m,n}(x, y) .

For the case y = 1 (zero magnetic field), this is known explicitly, thanks to the mathematical tour de force (using very abstract representation theory!) of 1968 Chemistry Nobelist, physicist Lars Onsager. For general y this is a wide-open, possibly impossible, problem. So all we can hope for, so far, is getting ‘close to infinity’. But, following Onsager, we can try and do the semi-inifite case, and get explicit expressions for

A_{m}(x, y) := \frac{1}{m} lim_{n→∞} n^{−1} log Z_{m,n}(x, y) .

It is well known [5] (and easy to see) that A_{m}(x, y) is the largest root of the characteristic equation of the transfer matrix T(m)

det(T(m) − λI) = 0 .

Hence, it is a solution of an algebaric equation (by an appropriate change of independent and dependent variables, it could be viewed as formal power series). It also follows that the logarithm (and of course, its derivatives with respect to both variables) are holonomic formal power series (or functions) in both variables.

8. How Onsager’s solution could have been conjectured via the C-finite ansatz

Onsager’s solution was an amazing tour de force using human ingenuity and lots of ad hoc tricks to find an ‘explicit’ expression for the largest eigenvalue of the transfer matrix for the k ×infinity case for any (symbolic) k. It turned out this (at most) algebraic function of degree 2^{k} can be expressed as some explicit product of k algebraic functions of degree 2. Using the methods in our previous paper [2], this fact could be deduced automatically, and we believe that with a little more effort, the actual Onsager factorization would be automatically derivable for k ≤ 9, from which it would be immediate to conjecture the Onsager factoriziation for symbolic (i.e. general k). Of course, this would not be a rigorous proof, in the mathematical sense, but definitely good enough for physicists.

9. How Onsager’s solution would be conjectured via the Holonomic ansatz

The exact value of the largest eigenvalue λ_{1}^{(m)} of the 2^{m} × 2^{m} transition matrix, is only a stepping stone for the free energy, the limit of log λ_{1}^{(m)}/m as m goes to infinity, viewed as a function of the continuous parameter x (right now are only considering the zero magnetic field case, i.e. y = 1). It follows from general theory, that for each m log λ_{1}^{(m)} is holonomic (aka D-finite) in the variable x, i.e. satisfies a linear differential equation with polynomial coefficients. Alas, as m gets larger and larger, these differential equations get more and more complicated, (to see them for 2 ≤ m ≤ 6) go to the webpage of this article. Suprisingly the case of interest m = ∞ satisfies a reasonale differential equation. As is standard in the theory (see [5], Ch. 6), it is convenient to perform the change of variable

v = \frac{x − x^{−1}}{x + x^{−1}} .

In terms of this the free energy satisfies the linear differential equation

(−13 v^{6} +70 v^{5} −251 v^{4} +16 v^{3} +69 v^{2} +2 v −5)\frac{d}{dv} f(v)−v(v −1)(55 v^{5} −309 v^{4} +328 v^{3} +220 v^{2}−

87 v + 1)\frac{d^{2}}{dv^{2}} f(v) − v^{2}(33 v^{4} − 172 v^{3} − 10 v^{2} + 116 v − 15)(v − 1)^{2} \frac{d^{3}}{dv^{3}} f(v)

−4 v^{3}v + 1)(v^{2} − 6 v + 1)(v − 1)^{3} \frac{d^{4}}{dv^{4}} f(v) = 0 .

References

[1] Johannes Blümlmlein, Mark Round, Carsten Schneider Refined Holonomic Summation Algorithms in Particle

    Physics, https://arxiv.org/abs/1706.03677

[2] Manuel Kauers and Doron Zeilberger Factorization of C-finite Sequences, in: “Advances in Computer Algebra In Honour of Sergei Abramov’s’ 70th Birthday”, WWCA 2016, Waterloo, Schneider, Carsten, Zima, Eugene (Eds.), Springer

[3] Per Hakan Lundow , Ising Data, http://www.theophys.kth.se/ phl/Data/,

[4] G. Nandhini, M. Vinoh Kumar, and M.V. Sangaranarayanan Estimation of Magnetization, SUsceptibility and Specific heat for the Two-Dimensional Ising model in a non-zero magnetic field,

    https://arxiv.org/abs/1006.0351

[5] Colin J. Thompson “Mathematical Statistical Mechanics”, Macmillan, NY, 1972.

[6] J.M. Yeomans “Statistical‘ Mechanics of Phase Transitions”, Oxford University Press, Oxford, 1992.

Manuel Kauers, Institute of Algebra Computation, J. Kepler University Linz, Austria

E-mail address: manuel dot kauers at jku dot at

Doron Zeilberger, Mathematics Department, Rutgers University (New Brunswick), Piscataway, NJ, USA

E-mail address: DoronZeil at gmail dot com