\documentclass[a4paper,draft]{amsart}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{bbm}
\usepackage{a4wide}

\parindent=0pt
\parskip\smallskipamount

\let\set\mathbbm
\def\<#1>{\langle#1\rangle}
\let\ideal\unlhd
\def\i{\mathrm{i}}
\def\e{\mathrm{e}}
\def\O{\mathrm{O}}
\def\Z{\mathcal{Z}}
\def\M{\mathcal{M}}
\def\F{\mathcal{F}}
\newtheorem{thm}{Theorem}
\newtheorem{prop}{Proposition}
\newtheorem{defn}{Definition}
\newtheorem{conj}{Conjecture}

\newcommand\todo[1][.]{\edef\tmpa{.}\edef\tmpb{#1}%
  \ifx\tmpa\tmpb
    \typeout{To Be on page \thepage}\fbox{\bf To Be}
  \else
    \typeout{To Be on page \thepage: #1}\fbox{{\bf To Be:} #1}
  \fi
}

\begin{document}

 \author[Manuel Kauers and Doron Zeilberger]{Manuel Kauers\,$^1$ and Doron Zeilberger}

 \address{Manuel Kauers, Institute of Algebra  Computation, J. Kepler University Linz, Austria}
 \email{manuel dot kauers at jku dot at}

 \address{Doron Zeilberger, Mathematics Department, Rutgers University (New Brunswick), Piscataway, NJ, USA}
 \email{DoronZeil at gmail dot com}

 \thanks{$^1$ Partially supported by the Austrian FWF grant Y464-N18.}

 \date{ Feb. 15, 2018}
 \title[Symbolic Computation to the Aid of Statistical Mechanics]
{Symbolic Computation to the Aid of Statistical Mechanics}

  \begin{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.

  \end{abstract}

 \maketitle

\section{Prologue}

Until about thirty years ago, physics consisted of {\bf two} parts: theoretical and experimental.
Then a new kid came to the block, {\it 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 {\it 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 \cite{brs}. Here we hope to illustrate the usefullness of symbolic computation to statistical mechanics, by focusing,
as a {\it case study}, on the celebrated {\it Ising model} (e.g. \cite{cjt72}).

\section{A crash course in statistical mechanics}

The general scenario of {\it statistical mechanics} consists of very many possible {\it configurations}, each with a
certain {\it 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  {\bf proportional} to
$$
e^{\frac{-E(C)}{kT}} \quad ,
$$
where $k$ is {\bf Boltzmann's constant}, and $T$ is the {\bf absolute temperature}.

It follows that the higher the energy, the less likely it is to show up, the 
decay being {\it exponential}.

In order to get {\it macroscopic} (thermodynamic)
(average) information, we form the famous {\bf Gibbs Partition function}
$$
\Z = \sum_{C \in configurations} e^{\frac{-E(C)}{kT}} \quad ,
$$
and the `average' energy (that physicists call {\it internal energy}) is obtained by
taking the {\it 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 {\it toy models} that `simplify' nature, and try to preserve
the {\it qualitative} features of the real world. Unfortunately, even these {\it toy} models
are extremely difficult, and with rare exceptions, intractable (or at least wide open), and
people resort to {\it approximations}, using {\bf numerical} methods, sophisticated brilliant
(but non-rigorous) renormalization theory calculations, and most often, {\it simulations} (aka ``Monte Carlo methods'').

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

\section{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 \times n$ matrices whose entries are $-1$ or $1$ (called `spin down' and `spin up', respectively).

Let
$$
x=e^{\frac{J}{kT}} \quad , \quad y=e^{\frac{H}{kT}}  \quad
$$
where $J$ is the so-called {\it coupling constant} (and hence a (physical) {\bf constant})
and $H$ is another input variable (in addition to $T$) called the (strength of the) {\it external magnetic} field.

The {\it energy} associated to a specific matrix 
$M \in \M(m,n)$ is defined as follows
$$
E(M):=-J \left ( \sum_{i=1}^{m} \sum_{j=1}^{n}  \left ( M_{i,j} M_{i,j+1} + M_{i,j} M_{i,j+1} \right )  \right )
-H \left ( \sum_{i=1}^{m} \sum_{j=1}^{n}  M_{i,j} \right ).
$$
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 {\it weight} of a matrix is
$$
weight(M):=
exp ( -\frac{E(M)}{kT} )=
x^{  \sum_{i=1}^{m} \sum_{j=1}^{n} \left ( M_{i,j} M_{i,j+1} + M_{i,j} M_{i,j+1} \right ) } \cdot
y^{  \sum_{i=1}^{m} \sum_{j=1}^{n}  M_{i,j}  } \quad .
$$

Our first (modest!) goal is just to compute, for as large $m$ and $n$ that our computers and human ingenuity
would allow, the (Gibbs) {\it Partition function} (i.e. {\it weight enumerator} of the
set $\M(m,n)$).
$$
\Z_{m,n}(x,y):=\sum_{M \in \M(m,n)} weight(M) \quad ,
$$
a certain {\it Laurent polynomial} in the variables $x$ (of degree $2mn$ and low-degree $ \geq -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=\infty$).

\section{Transfer matrices}

A {\it direct} \cite{nks10} 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 \cite{cjt72,jy},
called the {\it 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, \dots , 2b_{m-1}-1 )$, where $i=\sum_{j=0}^{m-1} b_j 2^j$ is the binary representation
of $i$.

Define the $2^m \times 2^m$ transfer matrix
$$
\mathbf{T}(m):=(T_{ij})_{0 \leq i,j  \leq 2^{m} -1} \quad,
$$
where the typical entry $T_{ij}$ is as follows.
$$
T_{ij}:=x^{\sum_{l=1}^m v(i)_l v(j)_l + \sum_{l=1}^{m-1} v(j)_l v(j)_{l+1} + v(j)_m v(j)_1} \cdot
y^{\sum_{l=1}^m v(j)_l} \quad .
$$

Using the transfer matrix, the desired polynomial 
$\Z_{m,n}(x,y)$ can be expressed very succinctly as
$$
\Z_{m,n}(x,y)= \mathrm{Tr} \,\, [ \, \mathbf{T}(m)^n \, ] \quad .
$$

It is well-known (and easy to see!) that {\it 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 {\it numeric} matrix.
Here we have a {\it symbolic} matrix, and a naive approach, (once $m$ gets larger)  would  be unfeasible.
But with the method of {\it homomorphic images}, we can go pretty far.


\section{The Functions that physicists care about}

For the finite (torodial) $m \times n$ lattice,
one is interested in the {\it Free energy},
$\F$
$$
\F= -kT \log \Z \quad ,
$$
{\it Internal Energy}, $U$\,
$$
U= kT^2 \frac{\partial}{\partial T} \log \Z \quad ,
$$
and {\it Specific heat}
$$
C= \frac{ \partial U}{\partial T} \quad .
$$
Physicists are also very interested in
the {\it Magnetization}, $M$
$$
M= -\frac{\partial \F}{\partial H} \quad ,
$$
and the {\it Isothermal susceptibilty}
$$
\chi_T = \frac{\partial M}{\partial H} \quad .
$$
At the ``end of the day'', people take the limit as
{\bf 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.

\section{ The Data for $\Z_{m,n}(x,y)$}

For $\Z_{n,n}(x,y)$  for $1 \leq n \leq 9$ see:
\begin{center}
\verb|http://www.math.rutgers.edu/~zeilberg/tokhniot/ising/Polys|
\end{center}



Very impressively, Per Hakan Lundow\cite{phl}, has already computed these for $n \leq 16$

For $\Z_{m,n}(x,y)$  for $1 \leq m \leq 5 \quad , \quad 1 \leq n \leq 11$ 
see:
\begin{center}
\verb|http://www.math.rutgers.edu/~zeilberg/tokhniot/ising/Polys1|
\end{center}





\section{ 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} \quad .
$$
But, physicists (and chemists)  that believe in the (fictional!) infinity, are {\it really} interested in
the {\it free energy} $\F(x,y)$
$$
\F:=
- kT \lim_{ m \rightarrow \infty}
\lim_{ n \rightarrow \infty} (mn)^{-1} \log \Z_{m,n}(x,y) \quad.
$$
For the case $y=1$ (zero magnetic field), this is known {\it explicitly}, thanks to
the {\bf mathematical} {\it tour de force} 
(using very abstract representation theory!) of 1968 {\bf Chemistry} Nobelist, {\bf 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 {\it explicit}
expressions for
$$
A_{m}(x,y) :=
\frac{1}{m} \lim_{ n \rightarrow \infty} n^{-1} 
\log \Z_{m,n}(x,y) \quad.
$$

It is well known \cite{cjt72} (and easy to see) that 
$A_m(x,y)$ is the largest root of
the {\it characteristic equation} of the 
transfer matrix $\mathbf{T}(m)$
$$
\det (\mathbf{T}(m)- \lambda I) = 0 \quad .
$$
Hence,  it is a solution of an {\it algebaric} equation (by an appropriate change of
independent and dependent variables, it could be viewed as {\it formal power series}).
It also follows that the {\it logarithm} (and of course, its derivatives with respect to both
variables) are {\it holonomic} formal power series (or  functions) in both variables.



\section{How Onsager's solution could have been conjectured via the C-finite ansatz}

Onsager's solution was an amazing {\it 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 \times infinity$ case for {\it 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  \cite{kz},
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 \leq 9$, from which it would be immediate
to {\bf conjecture} the Onsager factoriziation for {\it symbolic} (i.e. general $k$). Of course, this
would not be a rigorous proof, in the mathematical sense, but definitely good enough for physicists.

\section{How Onsager's solution would be  conjectured via the Holonomic ansatz}

The exact value of the largest eigenvalue $\lambda_1^{(m)}$ of the $2^m \times 2^m$ transition matrix,
is only a stepping stone for the free energy,
the limit of $\log \lambda_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 \lambda_1^{(m)}$ is {\it holonomic}
(aka $D$-finite) in the variable $x$, i.e. satisfies a linear differential equation with {\it polynomial} coefficients.

Alas, as $m$ gets larger and larger, these differential equations get more and more complicated, 
(to see them for $2 \leq m \leq 6$) go to the webpage of this article.

Suprisingly the case of interest $m=\infty$ satisfies a reasonale  differential equation.

As is standard in the theory (see \cite{cjt72}, Ch. 6), it is convenient to perform the change of variable
$$
v= \frac{x-x^{-1}}{x+x^{-1}} \quad .
$$
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}}{d{v}^{2}}}f ( v ) -{v}^{2} ( 33\,{v}^{4}-172\,{v}^{3}-10\,{v}^{2}+116\,v-15 )  ( v-1 ) ^{2}{\frac {d^{3}}{d{v}^{3}}}f ( v ) 
$$
$$
-4
\,{v}^{3}  v+1 )  ( {v}^{2}-6\,v+1 )  ( v-1 ) ^{3}{\frac {d^{4}}{d{v}^{4}}}f ( v) =0 \quad .
$$





\bibliographystyle{plain}
\bibliography{all}

\begin{thebibliography}{10}



\bibitem{brs} Johannes Bl\"umlmlein, Mark Round, Carsten Schneider
\emph{Refined Holonomic Summation Algorithms in Particle Physics},
{\tt https://arxiv.org/abs/1706.03677}


\bibitem{kz} Manuel Kauers and Doron Zeilberger
\emph{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


\bibitem{phl} Per Hakan Lundow ,
\emph{Ising Data, http://www.theophys.kth.se/~phl/Data/},

\bibitem{nks10} G. Nandhini, M. Vinoh Kumar, and M.V. Sangaranarayanan
\emph{Estimation of Magnetization, SUsceptibility and Specific heat for the Two-Dimensional Ising model in a non-zero magnetic field},
{\tt https://arxiv.org/abs/1006.0351} 


\bibitem{cjt72} Colin J. Thompson
\emph{``Mathematical Statistical Mechanics''},
Macmillan, NY, 1972.



\bibitem{jy} J.M. Yeomans
\emph{``Statistical` Mechanics of Phase Transitions''},
Oxford  University Press, Oxford, 1992.





\end{thebibliography}
 
\end{document}






