How did Ernst Ising originally solve the “Ising model?”
In 1920 Wilhelm Lenz proposed that ferromagnetism could be explained through the local interactions of spins in a material, and in1922 he gave this problem to his Ph.D. student Ernst Ising who ultimately solved it in one-dimension for his 1924 thesis. The hope was that the model could serve as an example of how local interactions could lead to non-local order, but this one-dimensional system didn’t exhibit the desired phase transition properties. Still, this first solution laid the foundation for higher-dimensional versions which validated Lenz’s original proposal (See Richard Folk’s “The Survival of Ernst Ising and the Struggle to Solve His Model” for a history of the early and later developments for the model).
Today, the Ising model is studied using a transfer-matrix formalism which greatly simplifies Ising’s original analysis. Still reviewing his original analysis is interesting due to its focus on the combinatorial representation of the partition function. It also makes one appreciate the slick-ness of the current formulation.
System Setup
Say that we have a system consisting of a linear chain of \(N\) spins each of which is either positive or negative. The system exists in a magnetic field \(H\) that contributes to the total energy with \(-\mu H\) for positive spins and \(+\mu H\) for negative spins. Adjacent spins that are different contribute an energy \(+\epsilon\) to the system; adjacent spins that are the same contribute no energy to the system.
We define a particular state by \(n_{\,+}\), \(n_{\,-}\), \(s\), and \(\delta\). The integers \(n_{\,+}\) and \(n_{\,-}\) represent, respectively, the total number of positive and negative spins in the system. The integer \(s\) is the number of negative-spin gaps between groups of positive spins, for systems that begins with a positive spin. For systems that begin and end with the same spin, we set \(\delta =0\), and for systems that begin and end with different spins we set \(\delta =1\). Take the following state as an example:
\begin{equation}
+ \overbrace{- \,-\, -} + + \overbrace{-\,-} + \overbrace{-\,-\,-} + -\,-\,- \quad (\text{State example})
\end{equation}
For this state, we have \(n_{\,+} = 6\) positive spins, \(n_{\,-} = 11\) negative spins, \(s = 3\) gaps between the groups of positive spins, and \(\delta = 1\) since the . We note that \(s+\delta\) is the number of “baskets” in which the negative spins are distributed, and \(s\) is the number of times we switch from a negative spin to a positive spin when moving from left to right.
Now let’s do some goal setting: Lenz’s ultimate objective in proposing this model was to explain ferromagnetism, the appearance of magnetization in systems that are not subject to a magnetic field. So we would like to compute the average magnetization for this collection of spins and see whether this magnetization is non-zero for zero-magnetic fields. To compute the average magnetization, we need to first compute the partition function since the average magnetization can be written in terms of the partition function’s derivative. Finally, to compute the partition function we must define the degeneracy factor for our space of states and the associated energy.
So here’s the sequence of steps we will hollow to reach the goal of computing the average magnetization
- Compute degeneracy factor for a state
- Ccompute energy for a state
- Compute partition function from energy and degeneracy
- Use partition function to compute magnetization
Let’s implement these steps in sequence!
Degeneracy Factor
For our system we define a state by the parameters \(n_{\,+}\), \(n_{\,-}\), \(s\), and \(\delta\). However, there are many such states for any particular choice of these parameters. So to properly perform the sum over states as it appears in the partition function, we must determine the number of different states are associated with the particular set of parameters \(n_{\,+}\), \(n_{\,-}\), \(s\), and \(\delta\).
First, say we only have \(n_{\,+}\) spins arranged in a line. How many ways can we select \(s\) locations to put gaps between these spins? Since there are \(n_{\,+}-1\) gaps between a list of \(n_{\,+}\) spins, we find that the number of ways to select \(s\) of them is
\begin{equation}
\binom{n_{\,+}-1}{s}. \qquad (1)
\end{equation}
Now, say that we have \(n_{\,-}\) identical objects and \(s + \delta\) different baskets in which we can place these objects. If each basket must contain at least one object, how many ways can we arrange the objects in baskets? (This is the number of ways to arrange the \(n_{\,-}\) negative spins amongst the \(s + \delta\) baskets available for these spins). This question is equivalent to asking how many ways can we select \(s + \delta -1\) locations to put gaps between \(n_{\,-}\) spins. Thus, given Eq.(1), we find that the answer to this second question is
\begin{equation}
\binom{n_{\,-}-1}{s+\delta-1}. \qquad (2)
\end{equation}
For states, that begin with positive spins, the product of Eq.(1) and Eq.(2) represents the number of different possible states that have the parameters \(n_{\,+}\), \(n_{\,-}\), \(s\), and \(\delta\), i.e., the degeneracy factor. Denoting this factor as \(\Omega_{\,+}\) we have
\begin{equation}
\Omega_{\,+} = \binom{n_{\,+}-1}{s}\binom{n_{\,-}-1}{s+\delta-1} \qquad (3)
\end{equation}
However, we can also have the reversed situation where we begin with negative spins rather than positive spins. In this case, all of our \(n_{\,+}\) change to \(n_{\,-}\) and vice versa. The degeneracy factor for this case, denoted \(\Omega_{\,-}\) is then
\begin{equation}
\Omega_{\,-} = \binom{n_{\,-}-1}{s}\binom{n_{\,+}-1}{s+\delta-1} \qquad (4)
\end{equation}
The net degeneracy factor is a combination of these two cases. Thus, we see that the number of different states with the parameters \(n_{\,+}\), \(n_{\,-}\), \(s\), and \(\delta\) is
\begin{equation}
\Omega(n_{\,+}, n_{\,-}, s, \delta ) = \binom{n_{\,+}-1}{s}\binom{n_{\,-}-1}{s+\delta-1}+\binom{n_{\,-}-1}{s}\binom{n_{\,+}-1}{s+\delta-1}. \qquad (5)
\label{eq:degen_tot}
\end{equation}
Energy
Next, in the task list is the system energy. Like the degeneracy, we want this energy to be a function of the parameters \(n_{\,+}\), \(n_{\,-}\), \(s\), and \(\delta\). Since there are \(s\) negative-spin gaps between positive spins, there are \(2s\) opposite-spin interactions from these gaps. There is an additional opposite-spin interaction for \(\delta=1\). Therefore, the interaction energy for the system is \(\epsilon (2s + \delta) \). The \(n_{\,+}\) positive spins and \(n_{\, -}\) negative spins, contribute an energy \(\mu H(n_{\,-}- n_{\,+})\) to the system. In all, then, the energy of the state is
\begin{equation}
E(n_{\,+}, n_{\,-}, s, \delta ) = \epsilon(2 s + \delta) + \mu H(n_{\,-}-n_{\,+}). \qquad (6)
\label{eq:en_tot}
\end{equation}
Partition Function
With Eq.(5) and Eq.(6), we can now write down the partition function for this system
\begin{align}
Z_N(\gamma, \alpha) &= \sum_{s=0}^{\infty} \sum_{\delta =0}^{1}\sum_{n_{\,+}=0}^{\infty} \sum_{n_{\,-} = 0}^{\infty}\left[\binom{n_{\,+}-1}{s}\binom{n_{\,-}-1}{s+\delta-1}+ \binom{n_{\,-}-1}{s}\binom{n_{\,+}-1}{s+\delta-1}\right]\\[.5em]
& \hspace{2.5cm} \times \delta_{N, n_{\,+}+n_{\,-}}\,e^{- (2s + \delta) \gamma + (n_{\,+}- n_{\,-}) \alpha}, \qquad (7)
\label{eq:ising_part}
\end{align}
where we defined \(\gamma \equiv \beta \epsilon\) and \(\alpha \equiv \beta \mu H\), and the Kroenecker delta \( \delta_{N, n_{\,+}+n_{\,-}}\) ensures that we only include terms for which the total number of spins is \(N\).
Generating Functions
Eq.(7) is quite complicated, so we would like to reduce it to something simpler. We make progress towards this simplification by introducing a generating function. Starting from the definition
\begin{equation}
F(x; \gamma, \alpha) \equiv \sum_{N=0}^{\infty}Z_N(\gamma, \alpha) x^N, \qquad (8)
\label{eq:Fx_def0}
\end{equation}
we find that, in computing \(F\), summing over \(N\) with the Kroenecker delta factor leads to \(N\) being set to \(n_{\,+} + n_{\,-}\). We are then left with the expression
\begin{align}
F(x; \gamma, \alpha)&= \sum_{s=0}^{\infty} \sum_{\delta =0}^{1}\sum_{n_{\,+}=0}^{\infty} \sum_{n_{\,-} = 0}^{\infty}\left[\binom{n_{\,+}-1}{s}\binom{n_{\,-}-1}{s+\delta-1}+ \binom{n_{\,-}-1}{s}\binom{n_{\,+}-1}{s+\delta-1}\right]\\[.5em]
& \hspace{2.5cm} \times x^{n_{\,+}+n_{\,-}}e^{- (2s + \delta) \gamma + (n_{\,+}- n_{\,-}) \alpha}. \qquad (9)
\label{eq:Fx_sec}
\end{align}
The second term in the brackets of Eq.(9) switches the roles of \(n_{\,+}\) and \(n_{\,-}\). Given the exponential multiplying this term, we see that relative to the first term in the brackets this switch is equivalent to taking \(\alpha \to -\alpha\). Therefore, after rearranging some terms, we find \(F\) can be expressed as
\begin{align}
F(x; \gamma, \alpha)&= \sum_{s=0}^{\infty} e^{-2\gamma s} \sum_{\delta =0}^{1} e^{-\gamma \delta} \sum_{n_{\,+}=0}^{\infty}\binom{n_{\,+}-1}{s} (x e^{\alpha})^{n_{\,+}} \sum_{n_{\,-}=0}^{\infty} \binom{n_{\,-}-1}{s+\delta-1} (x e^{-\alpha})^{n_{\,-}}\\[.5em]
& \hspace{2.5cm} + (\alpha \to -\alpha). \qquad (10)
\end{align}
From the identity \(\sum_{k=0}^{\infty}\binom{\ell-1}{k} x^k = (x/(1-x))^{\ell+1}\), we find
\begin{align}
\sum_{n_{\,+}=0}^{\infty}\binom{n_{\,+}-1}{s} (x e^{\alpha})^{n_{\,+}}& \sum_{n_{\,-}=0}^{\infty} \binom{n_{\,-}-1}{s+\delta-1} (x e^{-\alpha})^{n_{\,-}}\\[.5em]
& = \left(\frac{x e^{\alpha}}{1 – x e^{\alpha}} \right)^{s+1} \left( \frac{x e^{-\alpha}}{1 – x e^{-\alpha}} \right)^{s+\delta}\\[.5em]
& = \left(\frac{x^2}{1 -2x \cosh\alpha + x^2}\right)^{s} \left( \frac{x e^{-\alpha}}{1 – x e^{-\alpha}} \right)^{\delta}\left( \frac{x e^{\alpha}}{1 – x e^{\alpha}} \right). \qquad (11)
\end{align}
So, from the summation over \(\delta\), we have
\begin{align}
\left( \frac{x e^{\alpha}}{1 – x e^{\alpha}} \right)\sum_{\delta =0}^{1} \left( \frac{x e^{-\alpha-\gamma}}{1 – x e^{-\alpha}} \right)^{\delta} & = \left( \frac{x e^{\alpha}}{1 – x e^{\alpha}} \right)\frac{1- x e^{-\alpha}(1 – e^{-\gamma})}{1 – x e^{-\alpha}}\\[.5em]
& = \frac{x\big(e^{\alpha} – x(1-e^{-\gamma})\big)}{1- 2x \cosh \alpha + x^2} \qquad (12)
\end{align}
And the summation over \(s\) yields
\begin{equation}
\sum_{s=0}^{\infty} \left(\frac{x^2e^{-2\gamma}}{1 -2x \cosh\alpha + x^2}\right)^{s} = \frac{1- 2x \cosh \alpha + x^2}{1- 2x \cosh \alpha + (1- e^{-2\gamma})x^2}. \qquad (13)
\end{equation}
Therefore, \(F\) reduces to
\begin{align}
F(x; \gamma, \alpha) & = \frac{1- 2x \cosh \alpha + x^2}{1- 2x \cosh \alpha + (1- e^{-2\gamma})x^2} \cdot\frac{x\big(e^{\alpha} – x(1-e^{-\gamma})\big)}{1- 2x \cosh \alpha + x^2} + (\alpha \to\, – \alpha)\\[.5em]
& = \frac{x\big(e^{\alpha} – x(1-e^{-\gamma})\big)}{1- 2x \cosh \alpha + (1- e^{-2\gamma})x^2}+ (\alpha \to \,- \alpha), \qquad (14)
\end{align}
which upon simplification becomes
\begin{equation}
F(x; \gamma, \alpha) = \frac{2x\big(\cosh\alpha – x(1-e^{-\gamma}) \big)}{1 – 2x \cosh\alpha + x^2(1-e^{-2\gamma})}. \qquad (15)
\label{eq:Fx_res0}
\end{equation}
As a check of this result we can set \(\gamma = 0\) (corresponding to no interactions between spin) in Eq.(15), to obtain
\begin{equation}
F(x; \gamma= 0, \alpha) = \frac{2x\cosh\alpha}{1 – 2x \cosh\alpha } = \sum_{N=0}^{\infty} \left(2\cosh\alpha\right)^{N+1} x^{N+1}.
\end{equation}
which, given the definition of the generating function, implies that \(Z_{N}(\gamma = 0, \alpha) = (2 \cosh\alpha)^{N}\). This partition function value is consistent with the known paramagnetic result.
Expanding Generating Function
With our simplified generating function expression Eq.(15), we can now expand it and extract the coefficient of \(x^N\) to obtain the partition function. To achieve this expansion, we will have to modify Eq.(15) a bit.
First, we note that through the quadratic formula, we have
\begin{equation}
1 – 2x \cosh\alpha + x^2(1-e^{-2\gamma}) = (1 – e^{-2\gamma}) (x_+- x)(x_- -x) \qquad (16)
\end{equation}
where
\begin{equation}
x_{\pm} \equiv \frac{\cosh\alpha \pm \sqrt{\sinh^2 \alpha +e^{-2\gamma}}}{1- e^{-2\gamma}}. \qquad (17)
\end{equation}
Also expanding the inverse product \(1/(x_+- x)(x_- -x)\) in \(x\), we obtain
\begin{align}
\frac{1}{(x_+- x)(x_- -x)} & = \frac{1}{x_+ – x_-} \left[ \frac{1}{x_+ – x} – \frac{1}{x_- – x}\right]\\[.5em]
& = \frac{1- e^{-2\gamma}}{2x\sqrt{\sinh^2 \alpha +e^{-2\gamma}}}\sum_{N=0}^{\infty}x^{N+1}\left[\left(\frac{1}{x_+}\right)^{N+1} -\left( \frac{1}{x_+}\right)^{N+1} \right]. \qquad (18)
\end{align}
Thus Eq.(15) becomes
\begin{equation}
F(x; \gamma, \alpha) = \frac{\cosh\alpha – x(1-e^{-\gamma})}{\sqrt{\sinh^2 \alpha +e^{-2\gamma}}}\sum_{N=0}^{\infty}x^{N+1}\left[\left(\frac{1}{x_+}\right)^{N+1} -\left( \frac{1}{x_-}\right)^{N+1} \right]. \qquad (19)
\label{eq:Fx_res2}
\end{equation}
Next, we define
\begin{equation}
\lambda_{\,\pm} \equiv \frac{1}{x_\mp} = \cosh \alpha \pm \sqrt{\sinh^2 \alpha +e^{-2\gamma}}, \qquad (20)
\end{equation}
\begin{equation}
A \equiv \frac{\cosh\alpha }{\sqrt{\sinh^2 \alpha +e^{-2\gamma}}}, \qquad B \equiv \frac{1-e^{-\gamma}}{\sqrt{\sinh^2 \alpha +e^{-2\gamma}}}, \qquad (21)
\end{equation}
and
\begin{equation}
\bar{f}(N, \gamma, \alpha) \equiv \lambda_{\,+}^{N} – \lambda_{\,-}^{N}. \qquad (22)
\label{eq:fN_def}
\end{equation}
With these definitions Eq.(19) becomes
\begin{equation}
F(x; \gamma, \alpha) = (A – Bx ) \sum_{N=0}^{\infty} \bar{f}(N+1, \gamma, \alpha)x^{N+1}. \qquad (23)
\end{equation}
And with some additional rearranging, we find
\begin{equation}
F(x; \gamma, \alpha) = A \bar{f}(1, \gamma, \alpha)x + \sum_{k=1}^{\infty}x^{N+1} \Big[A\bar{f}(N+1, \gamma, \alpha) – B\bar{f}(N, \gamma, \alpha)\Big]. \qquad (24)
\end{equation}
From the original generating function expression Eq.(8), we can conclude that the partition function is (for \(N> 0\))
\begin{equation}
Z_{N+1}(\gamma, \alpha) = A\bar{f}(N+1, \gamma, \alpha) – B\bar{f}(N, \gamma, \alpha). \qquad (25)
\end{equation}
Substituting Eq.(22) into this expression leads to the final result
\begin{align}
Z_{N+1}(\gamma, \alpha) = p_{\,+}\lambda_{\,+}^{N}\, – \,p_{\,-}\lambda_{\,-}^{N}, \qquad (26)
\label{eq:ising_part0}
\end{align}
where
\begin{align}
p_{\pm} & \equiv A \lambda_{\pm}- B \\[.5em]
& = \frac{\cosh\alpha }{\sqrt{\sinh^2 \alpha +e^{-2\gamma}}}\left(\cosh \alpha \pm \sqrt{\sinh^2 \alpha +e^{-2\gamma}}\right) – \frac{1-e^{-\gamma}}{\sqrt{\sinh^2 \alpha +e^{-2\gamma}}}\\[.5em]
& = \frac{\sin^2\alpha + e^{-\gamma}}{\sqrt{\sinh^2 \alpha +e^{-2\gamma}}}\pm \cosh\alpha. \qquad (27)
\end{align}
Magnetization
With the partition function Eq.(26), we are now at last ready to compute the magnetization. We will make an approximation to perform this calculation. For a system in which \(N\gg1\), we find \(\lambda_+^N \gg \lambda_-^N\), and so the partition function Eq.(26) can be approximated as
\begin{equation}
Z_{N}(\gamma, \alpha) \simeq p_{\,+} \lambda_{\,+}^N. \qquad (28)
\end{equation}
Computing the mean magnetization in this same \(N\gg1\) limit, we have
\begin{align}
m & = \mu \frac{\partial}{\partial \alpha} \ln Z_N(\gamma, \alpha) \\[.5em]
& \simeq \mu N\frac{\partial}{\partial \alpha} \ln \left(\cosh \alpha \pm \sqrt{\sinh^2 \alpha +e^{-2\gamma}}\right)\\[.5em]
& = \frac{\mu N \sinh \alpha}{\sqrt{\sinh^2 \alpha +e^{-2\gamma}}},
\label{eq:ising_m} \qquad (29)
\end{align}
where we dropped the \(\ln p_+\) term in the second line because it is subleading in the \(N\gg1\) limit. Eq.(29) is the standard expression for the magnetization and is the result first found by Ising in his thesis1.
(Figure 1: English translation of result from Ising’s Thesis. Link in footnotes)
As noted by Ising, the important thing to note is that in the \(\alpha \to 0\) (i.e., zero magnetic field) limit, the magnetization \(m\) Eq.(29) goes to zero, implying this system does not exist ferromagnetism.
Final Remarks
We just went through Ernst Ising’s original calculation of the partition function and magnetization of a system of one-dimensional spins. The partition function calculation involved careful combinatorial arguments, generating functions, and expansions and was, in general, much more complicated than the modern derivation using transfer matrices. To compare, consider a sketch of the modern derivation. We represent the \(N\) spin system by the vector \(s_1, s_2, \ldots, s_N\) where \(s_i \in \{-1, +1\}\). The partition function for the system is then
\begin{align}
Z_N(\gamma, \alpha) & = \sum_{s_1, s_2, \ldots, s_N} \exp\left(-\frac{\gamma}{2}\sum_{i=1}^{N-1}(1-s_{i}s_{i+1}) + \alpha \sum_{i=1}^N s_i\right) \\
& = \sum_{s_1, s_2, \ldots, s_N} V_{s_1, \,s_2} V_{s_2, \,s_3} \cdots V_{s_{N-1}, \,s_N}\\
&=\sum_{s_1, s_N} (V^{N-1})_{s_1, \, s_N} \qquad (30)
\end{align}
where we defined the transfer matrix \(V\) with elements \( V_{s_i, \, s_j} = \exp\left(-\frac{\gamma}{2}(1-s_{i}s_{j}) + \frac{\alpha}{2} (s_i + s_j)\right)\) and with matrix representation
\begin{equation}
V = \left(\begin{array}{c c}V_{+1, +1} &V_{+1, -1} \\ V_{-1, +1} & V_{-1, -1} \end{array}\right)= \left(\begin{array}{c c}e^{\alpha} & e^{-\gamma} \\ e^{-\gamma} & e^{-\alpha} \end{array}\right). \qquad (31)
\end{equation}
Next, we diagonalize \(V\) with the change of basis formula \(V = O \Lambda O^T\) for which \(\Lambda = \text{diag}(\lambda_{\,+}, \lambda_{\,-})\), where \(\lambda_{\, \pm}\) are the eigenvalues of \(V\) and \(O = (\textbf{v}_{+}, \textbf{v}_{-})\) is the matrix of normalized eigenvectors of \(V\).
\begin{align}
Z_N(\gamma, \alpha) &=\sum_{s_1, s_N}(O \Lambda^{N-1} O^T)_{s_1, \, s_N}, \\
&= c_{\,+} \lambda_{\, +}^{N-1}+ c_{\,-} \lambda_{\, +}^{N-1},
\end{align}
where \(c_{\,\pm}\) are coefficients determined from the elements of change of basis matrix \(O\). Thus, we are able to reproduce the partition function result Eq.(26) in much fewer steps than in Ising’s original derivation. This drastic reduction in derivation complexity is a good lesson in the importance of seeking better mathematical representations of problems.
Footnotes
- https://www.hs-augsburg.de/~harsch/anglica/Chronology/20thC/Ising/isi_fm00.html âŠī¸