Spin Glass Models 3: Ising Model - Theory
In this blog post we will introduce another model of spin glasses: the Ising model. We relax some of the simplifications of previous models to create something that more accurately captures the structure of real spin-glasses. We will look at this model mathematically to gain insight into its properties.
- Introducing the Ising Model
- 1D Ising Model
- 2d Ising Model
- Infinte Dimension Ising Model
- n-d Ising Model ($n\geq3$)
- Conclusion
This is the third blog post in a series - you can find the previous blog post here
In this blog post we are going to look at another spin-glass model. In the previous post we looked at the Sherrington-Kirkpatrick model which allowed us to study the spin-glass analytically. However there is a certain lack of "realism" in this model - the infinite interaction range means that the Sherrington-Kirkpatrick (in a sense) does not occupy any space, there is no concept of dimension nor geometry. While providing interesting mathematical results, some of which appear to be universal to spin-glass systems, we would like to look at a different model that captures more of the real world system and perhaps might uncover new behaviours.
Introducing the Ising Model
One model that we can look at is the Edwards-Anderson Ising model. We have actually looked at these models achronologically: Edwards-Anderson proposed the Ising model before Sherrington-Kirkpatrick proposed theirs. In fact Sherrington-Kirkpatrick developed their model partly due to the difficulty in deal with the Ising model. I have presented the models in this order in this series of blog posts because it feels more natural to look at models in increasing complexity order rather than presenting a historical account.
The main difference between Ising and Sherrington-Kirkpatrick is the extent over which the interactions can occur, instead of an infinite range the Ising model only allows "nearest-neighbour" interactions. This relaxation means that there is a concept of dimension and geometry to an Ising spin-glass. For example we could think of spins oriented on a line segment (a finite 1d Ising model), a square lattice (a 2d Ising model) or something more exotic like spins oriented on the nodes of a 32 dimensional hexagonal lattice. Through careful use of limits we could also consider infinite-dimensional Ising models too. The Hamiltonian follows the form one would expect: $$ H = - \sum_{\lt x,y \gt} J_{xy} \sigma_x \sigma_y - h \sum_x \sigma_x $$ Where we use the notation $\lt x,y \gt$ to denote the sum occuring over neighbour pairs.
As before the selection of $J_{xy}$ is somewhat arbitrary, however unlike the Sherrington-Kirkpatrick we do not have to scale this with lattice size to retain meaningful "average energy" or "average magnetism" measurements when we increase the size of the system. If we take $J_{xy} = J$ for some fixed $J$ we get the Ising model of Ferromagnetism, if we allow $J_{xy}$ to vary according to some distribution we get the Edwards-Anderson Ising Model. In this blog we will use the term "Ising model" to refer to either situation, the mathematics of moving from the fixed $J$ case usually involves integrating against the density function, this can lead to more complicated formulae so in this blog we will mainly focus on the ferromagnetic case unless otherwise stated.
As before spins can only take values $\pm 1$ - in some literature this is referred to as "Ising spins" (e.g. you can read references to "Sherrington-Kirkpatrick with Ising spins" - this means an infinite interaction range with binary spins).
A pictorial representation of a 3d square lattice Ising spin-glass can be seen below:
For the rest of this article we will limit ourselves to talking about "square" lattices since this is where the majority of research is focussed. Using different lattice types won't change the story too much.
In studying these models it is also worth noting what happens on the "boundary" of the lattice: for finite size systems this can become an issue. In this blog post we will skirt over the issue, where necessary we will assume periodic boundary conditions (i.e. there is a toroidal type geometry to the system) so each "edge" of the lattice is connected to an opposing one, so in a sense there are no boundary conditions to specify.
1D Ising Model
We will now consider a 1 dimensional Ising model where spins are located on a line segment. This simplication means that finding the ground state (minimal energy spin configuration) is trivial: we pick a site in the middle of the line segment (the origin) we will pick a spin ($\pm 1$) at random. We then propagate out from the origin in the left and right direction, we pick the spin for the next site according to the interaction strength for example: if $\sigma_0 = +1$ and $J_{0,1} > 0$ then $\sigma_1 = +1$ and so on. We can see that (as with all spin glass systems) there is a symmetry: if we flip all the spins direction we end up with an equivalent energy - we call these "pairs" of configurations. It doesn't take much convincing to realise that the pair of configurations we constructed is a unique minima for the system. To see this imagine we have only one pair of spins that opposes the direction indicated by the interaction strength. For a large enough system we can eventually find another interaction strength that is greater in magnitude (by definition having a spin pair as indicated by the sign of the interaction strength). By simply flipping the relative orientation of the pairs we will end up with a configuration of lower energy - which is a contradiction.
We will now look to solve the 1d Ising model analytically. For simplicity we will assume no external magnetic field, this just makes the formulae look "prettier". The crux of understanding any spin glass system is finding the Gibbs distribution, as before: $$P(\sigma) = \frac{exp(-\beta H_{\sigma})}{Z_{\beta}}$$ Where: $$H = - \sum_{\lt x,y \gt} J_{xy} \sigma_x \sigma_y$$ And $\beta$ is the inverse temperature. We can write an expression for the partition function as: $$ Z_{\beta} = \sum_{\sigma} exp(-\beta H_{\sigma}) = \sum_{\sigma} exp(-\beta \sum_{\lt x,y \gt} J_{xy} \sigma_x \sigma_y) $$ Suppose we look at a line segment with $N$ spins, we can then write this as: $$ Z_{\beta} =\sum_{\sigma_1, ..., \sigma_N} exp(-\beta (J_{1,2}\sigma_1 \sigma_2 + J_{2,3}\sigma_2 \sigma_3 + ... + J_{N-1,N}\sigma_{N-1}\sigma_N))$$ We notice that we can factorise this summation as: $$ Z_{\beta} =\sum_{\sigma_1, ..., \sigma_{N-1}} exp(-\beta (J_{1,2}\sigma_1 \sigma_2 + ... + J_{N-2,N-1}\sigma_{N-2}\sigma_{N-1}))\sum_{\sigma_N}exp(-\beta J_{N-1,N} \sigma_{N-1} \sigma_N) $$ The internal sum can then be evaluated: $$ \sum_{\sigma_N}exp(-\beta J_{N-1,N} \sigma_{N-1} \sigma_N) = exp(\beta J_{N-1,N} \sigma_{N-1} ) + exp(-\beta J_{N-1,N} \sigma_{N-1} ) = 2cosh(\beta J_{N-1,N} \sigma_{N-1}) = 2cosh(\beta J_{N-1,N}) $$ The last equality making use of the fact $\sigma_{N-1} = \pm 1$ and that $cosh(x) = cosh(-x)$. By repeating this process we can express the partition function as: $$ Z_{\beta} = 2 \prod_{i=1}^{N-1} 2 cosh(\beta J_{i,i+1}) $$ We can evaluate this exactly. If we are sampling the interaction strengths according to some distribution we can find expected values (or other statistical properties) by integrating against the density function as usual, since this adds to notational complexity we will assume that the interaction strengths are fixed for now. We can then calculate thermal properties using the equations:
\begin{align} F &= - \frac{1}{\beta} ln Z_{\beta} = - T ln 2 - T \sum_{i=1}^{N-1} ln (2 cosh(\beta J_{i,i+1})) \\ U &= - \frac{\partial}{\partial \beta} ln Z_{\beta} = - \sum_{i=1}^{N-1} J_{i,i+1} tanh(\beta J_{i,i+1}) \\ C &= \frac{\partial U}{\partial T} = \sum_{i=1}^{N-1} (\beta J_{i,i+1})^2 sech^2(\beta J_{i,i+1}) \\ S &= \frac{U - F}{T} = ln2 + \sum_{i=1}^{N-1} \left( -\beta J_{i,i+1} tanh(\beta J_{i,i+1}) + ln(2cosh(\beta J_{i,i+1}) \right) \end{align}Where:
F - is the Helmholtz-Free Energy
U - is the thermodynamic energy (ensemble average)
C - is the heat capacity
S - is the entropy
We can find the expected value of a given instantiation of interactions through an integral such as: $$ \mathbb{E}(U) = - \sum_{i=1}^{N-1} \int J_{i,i+1} tanh(\beta J_{i,i+1}) Q(J_{i,i+1}) dJ_{i,i+1} $$ With $Q(J)$ being the density function of the distribution and the integral occuring over its support. Similar formulae exist for the other thermodynamic variables and you can calculate other statistics (e.g. variance) in the usual way.
We can plot the values of these variables for a given set of interaction weights:
# This code creates a 4 figure plot showing how
# Thermodynamic variables change with temperature
# For a 1d Ising Model with Gaussian interactions
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# Fix seed
np.random.seed(123)
# Fix number of points
N = 100
# Instantiate interaction strength
J = np.random.normal(0,1,size=N-1)
# Set temperature ranges
T_min = 1e-4
T_max = 5
T_steps = 1000
T = np.arange(T_steps)/T_steps *(T_max - T_min) + T_min
beta = 1 / T
# Set up holders for variables
F = np.zeros(T_steps)
U = np.zeros(T_steps)
C = np.zeros(T_steps)
S = np.zeros(T_steps)
# Loop over T_steps and calculate at each step
for i in range(T_steps):
F[i] = - T[i] * np.log(2) - T[i] * (np.log(2*np.cosh(beta[i]*J))).sum()
U[i] = - (J * np.tanh(J * beta[i])).sum()
C[i] = ((beta[i] * J)**2 * (np.cosh(beta[i] * J))**-2).sum()
S[i] = (U[i] - F[i]) / T[i]
# Divide by number of points to give a scale invariant measure
F = F / N
U = U / N
C = C / N
S = S / N
# Create plots
fig, axs = plt.subplots(2, 2, figsize=(10,10), gridspec_kw={'hspace': 0.25, 'wspace': 0.25})
axs[0, 0].plot(T, F)
axs[0, 0].set_title("Free Energy")
axs[0, 0].set(xlabel="T", ylabel="F / N")
axs[0, 1].plot(T, U)
axs[0, 1].set_title("Average Energy")
axs[0, 1].set(xlabel="T", ylabel="U / N")
axs[1, 0].plot(T, C)
axs[1, 0].set_title("Heat Capacity")
axs[1, 0].set(xlabel="T", ylabel="C / N")
axs[1, 1].plot(T, S)
axs[1, 1].set_title("Entropy")
axs[1, 1].set(xlabel="T", ylabel="S / N")
plt.show()
From these plots we can see there is no phase transition taking place - this is true for all 1d Ising models.
We can also calculate the correlation function for the system. By following a similar line of logic as before we find: $$ \langle \sigma_n \sigma_{n+r} \rangle = \prod_{i=0}^{r} tanh(\beta J_{n+i, n+i+1})$$
We can plot this as a function of $r$ starting at the first position:
# This code creates a plot displaying
# The correlation function
# For a 1d Ising Model with Gaussian interactions
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# Fix seed
np.random.seed(123)
# Fix number of points
N = 100
# Instantiate interaction strength
J = np.random.normal(0,1,size=N-1)
# Create holders for variables
r = np.arange(N)
corr_array = np.zeros(N)
# Fix beta
beta = 1
# Calculate correlations
tanh_array = np.tanh(beta * J)
for i in range(N):
corr_array[i] = tanh_array[:i].prod()
plt.plot(r[0:20], corr_array[0:20])
plt.title("Correlation Function")
plt.ylabel(r"$\langle \sigma_{1} \sigma_{r} \rangle$")
plt.xlabel("r")
plt.xticks(np.arange(11)*2)
plt.ylim((-1,1))
plt.xlim((0,20))
plt.show()
As we can see there is a specific "structure" to the correlation function given an instantiation - this is not particularly surprising. If we picked an alternate starting state (e.g. the 2nd site) this graph can look totally different, the decay in absolute value will be similar however. We also notice that the correlation decays to zero fairly rapidly suggesting there isn't a long range structure to the model.
2d Ising Model
We now turn our attention to the 2d Ising model. The mathematics here will, understandably, get more complicated. We will require a bit more sophistication to solve this system. The solution was originally published by Lars Onsager in 1944, with alternative proofs and derivations published later. The derivation itself is fairly involved and would take a blog post (at least) by itself to cover - I may come back to this at a later date. For now I will simply present the result for the free energy (F) in absence of a magnetic field ($h=0$): $$ - \frac{\beta}{N} F = ln2 + \frac{1}{8\pi^2} \int^{2\pi}_{0} \int^{2\pi}_{0} ln \left[ cosh(2 \beta J_H) cosh(2 \beta J_V) - sinh(2 \beta J_H)cos(\theta_1) - sinh(2 \beta J_V)cos(\theta_2) \right]d\theta_1 d\theta_2 $$ Where instead of $J_{xy}$ being sampled from a gaussian distribution we assume fixed interaction strengths $J_H$ and $J_V$ in the horizontal and vertical directions.
For simplicity we will take: $J_H = J_V = J$ and we can derive the thermodynamic properties (per site - e.g. $f=\frac{F}{N}$) of this system as: \begin{align} f &= \frac{-ln2}{2\beta} - \frac{1}{2\pi\beta} I_0(\beta, J) \\ u &= -J coth(2\beta J) \left[ 1 + \frac{2}{\pi} \left( 2tanh^2(2\beta J) -1 \right) I_1(\beta, J) \right] \\ c &= 2J\beta^2 \left[U csch(2\beta J)sech(2\beta J) + \frac{8J}{\pi} sech^2(2\beta J) I_1(\beta, J) - \frac{2 \beta J}{\pi}(cosh(4\beta J)-3)^2 sech^6(2\beta J) I_2(\beta, J) \right] \\ s &= \frac{U - F}{T} \end{align}
For convenience I created 3 new functions $I_0(\beta,J), I_1(\beta, J)$ and $I_2(\beta, J)$ to make the equations a little shorter. These functions are defined as: \begin{align} I_0(\beta, J) &= \int^\pi_0 ln \left[ cosh^2(2 \beta J) + sinh^2(2 \beta J) \sqrt{1 + csch^4(2 \beta J) - 2 csch^2(2 \beta J) cos(2\theta)} \right] d\theta \\ I_1(\beta, J) &= \int^{\frac{\pi}{2}}_0 \left[1 - 4csch^2(2\beta J)( 1 + csch^2(2\beta J))^{-2} sin^2(\theta) \right]^{-\frac{1}{2}} d\theta \\ I_2(\beta, J) &= \int^{\frac{\pi}{2}}_0 sin^2(\theta)\left[1 - 4csch^2(2\beta J)( 1 + csch^2(2\beta J))^{-2} sin^2(\theta) \right]^{-\frac{3}{2}} d\theta \end{align}
As before we can produce plots of these:
# This code creates a 4 figure plot showing how
# Thermodynamic variables change with temperature
# For a 2d Ferromagnetic Ising Model with fixed interaction
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
%matplotlib inline
# Fix seed
np.random.seed(123)
# Instantiate interaction strength
J = 1
# Set temperature ranges
T_min = 1e-4
T_max = 5
T_steps = 1000
T = np.arange(T_steps)/T_steps *(T_max - T_min) + T_min
beta = 1 / T
# Set up holders for variables
f = np.zeros(T_steps)
u = np.zeros(T_steps)
c = np.zeros(T_steps)
s = np.zeros(T_steps)
# Set up integrands for I0, I1, I2
def integrand0(x, b, j):
return np.log(np.cosh(2*b*j)**2 + np.sinh(2*b*j)**2 * np.sqrt(1 + np.sinh(2*b*j)**(-4) - 2*np.sinh(2*b*j)**(-2) * np.cos(2*x)))
def integrand1(x, b, j):
return (1 - 4*np.sinh(2*b*j)**(-2)*(1 + np.sinh(2*b*J)**(-2))**(-2)*np.sin(x)**2)**(-0.5)
def integrand2(x, b, j):
return np.sin(x)**2 * integrand1(x, b, j)**3
# Loop over T_steps and calculate at each step
for i in range(T_steps):
bt = beta[i]
I0 = quad(integrand0, 0, np.pi, args=(bt, J))
I1 = quad(integrand1, 0, np.pi/2, args=(bt, J))
I2 = quad(integrand2, 0, np.pi/2, args=(bt, J))
f[i] = - np.log(2) / (2 * bt) - I0[0] / (2 * np.pi * bt)
u[i] = -J*np.tanh(2*bt*J)**(-1) * ( 1 + (2/np.pi)*(2*np.tanh(2*bt*J)**2 -1)*I1[0])
c[i] = 2*bt**2*J*( u[i]*np.sinh(2*bt*J)**(-1)*np.cosh(2*bt*J)**(-1) + (8*J / np.pi)*np.cosh(2*bt*J)**(-2)*I1[0] - (2*bt*J/np.pi)*((np.cosh(4*bt*J)-3)**2)*np.cosh(2*bt*J)**(-6)*I2[0] )
s[i] = (u[i] - f[i])*bt
# Create plots
fig, axs = plt.subplots(2, 2, figsize=(10,10), gridspec_kw={'hspace': 0.25, 'wspace': 0.25})
axs[0, 0].plot(T, f)
axs[0, 0].set_title("Free Energy")
axs[0, 0].set(xlabel="T", ylabel="f")
axs[0, 1].plot(T, u)
axs[0, 1].set_title("Average Energy")
axs[0, 1].set(xlabel="T", ylabel="u")
axs[1, 0].plot(T, c)
axs[1, 0].set_title("Heat Capacity")
axs[1, 0].set(xlabel="T", ylabel="c")
axs[1, 1].plot(T, s)
axs[1, 1].set_title("Entropy")
axs[1, 1].set(xlabel="T", ylabel="s")
plt.show()
We find there is a critical temperature $T_c$ where the specific heat equation diverges. We can compute the value of this as satisfying: $$ sinh\left(\frac{2J_H}{T_c}\right)sinh\left(\frac{2J_V}{T_c}\right) = 1 $$ In the case where $J_H = J_V = J$ we have: $$T_c = \frac{2 J}{ln\left(1+\sqrt{2}\right)} $$ This is an example of a second order phase transition since the discontinuity only arises under a second derivative.
For temperatres under this critical temperature we have that the spontaneous magnetization can be calculated as: $$ m = \left[ 1 - csch^2\left(\frac{2J_H}{T}\right)csch^2\left(\frac{2J_V}{T}\right) \right]^{\frac{1}{8}} $$
We can plot this (for $J_H = J_V = J = 1$) as:
# This code plots the spotaneous magnetization of
# a 2d Ising model on a square-lattice
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# Set up temperature array
T_min = 1e-4
T_max = 5
T_steps = 1000
T = np.arange(T_steps)/T_steps *(T_max - T_min) + T_min
# Fix interaction strength constant
J = 1
m = np.power(np.maximum(1 - np.sinh(2*J/T)**-4,0), 1/8)
plt.plot(T,m)
plt.xlabel("T")
plt.ylabel("Magnetization")
plt.title("2d Ising Model Spontaneous Magnetization")
plt.xlim((0,5))
plt.show()
print("Critical Temp (Tc):", (2*J) / np.log(1 + np.sqrt(2)))
The correlation function is harder to compute, here we just present the correlation function between elements on the diagonal of the lattice:
$$ \langle \sigma_{0,0} \sigma_{N,N} \rangle = det \left[ A_N \right]$$ Where $A_N = (a_{n,m})_{n,m=1}^N$ is an $NxN$ matrix with entries: $$ a_{n,m} = \frac{1}{2 \pi} \int_0^{2\pi} e^{i(n-m)\theta} \sqrt{\frac{sinh^2(2\beta J) - e^{-i\theta}}{sinh^2(2\beta J) - e^{i\theta}} } d\theta $$
We can plot this as so:
# This code plots the diagnoal correlation function
# For the 2d square-lattice Ising model
# Lattice size = NxN
# At the critical temperature
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
%matplotlib inline
# Set up integrand function
imag = complex(0,1)
def integrand(x, n, m, beta, J):
return np.exp(imag*(n-m)*x)*np.sqrt((np.sinh(2*beta*J)**2-np.exp(-imag*x))/(np.sinh(2*beta*J)**2-np.exp(imag*x)))
# Set up constants
J = 1
beta = np.log(1 + np.sqrt(2)) / (2*J)
# Set up arrays
N_max = 20
N_array = np.arange(N_max+1)
correl_array = np.zeros(N_max+1)
for x in range(N_max + 1):
N = N_array[x]
A_N = np.zeros((N, N))
for n in range(N):
for m in range(N):
I = quad(integrand, 0, 2*np.pi, args=(n, m, beta, J))
A_N[n,m] = I[0]/(2*np.pi)
correl_array[x] = np.linalg.det(A_N)
# Plot correlations
plt.plot(N_array, correl_array)
plt.xlabel("N")
plt.ylabel(r"$\langle \sigma_{0,0} \sigma_{N,N} \rangle$")
plt.title("Correlation Function")
plt.xticks(np.arange(11)*2)
plt.ylim((0,1))
plt.xlim((0,20))
plt.show()
We can see that even though we have short range interactions (nearest neighbour) this gives rise to longer-range structure within the model. In fact we find for temperatures below the critical temperature there is an infinite correlation length, whereas for temperatures above the critical temperature the correlation length is finite.
Infinte Dimension Ising Model
Next we consider the situation where the Ising model exists in infinite dimensions. In this case each site has infinitely many neighbours, as such a mean-field approximation is valid and we have a model that is somewhat similar to the Sherrington-Kirkpatrick fully connected geometry. (Please excuse the lack of rigour here; one has to be careful in how limits are defined and it is a non-trivial exercise. For this blog post I don't want to go down into that sort of detail.)
If each site has infinitely many neighbours we only need to concern ourselves with the ratio of positive and negative spin neighbours. By mean field if we take the probability of a positive spin as $p$ then via the Gibbs distribution we have: $$ \frac{p}{1-p} = exp(2\beta H)$$ The average magnetization can then be calculated: $$ \mathbb{E}\left[M\right] = (1)p + (-1)(1-p) = 2p - 1 = \frac{1 - exp(2\beta H)}{1 + exp(2\beta H)} = tanh(2 \beta H)$$
By investigating this function we can gain insight into spontaenous magnetization. Other similar arguments can be invoked for the other themodynamic properties. It is possible to show, as with the Sherrington-Kirkpatrick, that a phase transition occurs.
n-d Ising Model ($n\geq3$)
Finally we look at the case where the dimension of the model is finite but strictly greater than 2. In this situation things get much trickier, in fact there are not many defined mathematical results in these systems and this is the subject of current research. As such I will just briefly outline some of the approaches that have been suggested to study these systems and their properties (presented in chronological order from when they were proposed):
- Replica-Symmetry Breaking - From our previous post on Sherrington-Kirkpatrick models we briefly looked at this sort of method. They were an obvious first choice in trying to deal with short range interactions. It has been shown that the "standard" techniques are insufficient but there have been extensions proposed that have shown some promise. Like in the infinite range model it suggests states have an ultrametric structure and uncountably many pure states.
- Droplet Scaling - The first such argument being presented by Rudolf Peierls. The main idea is to consider the arisal of "loops" or "islands" of spins (clusters of atoms all with the same spin being enclosed by some boundary). We then aim to count the number of such loops, often in high or low temperature ranges via approximation. This leads to only 2 pure states.
- Chaotic Pairs - Has properties somewhat similar to a combination of the preceeding 2 methods, like replica symmetry breaking there are infinitely many thermodynamic states however the relationship between them is much simpler and has simple scaling behaviour - much like droplet scaling.
- TNT - This interpretation does not itself specify the number of pure states, however it has been argued that it most naturally exists with 2 pure states - much like droplet scaling. However it has scaling properties more similar to that of replica symmetry breaking.
For completeness: a pure state of a system is a a state that cannot be expressed as a convex combination of other states.
As far as I am aware there is currently no conesensus on which (if any) of the options presented is the correct interpretation.
One of the main questions that we would like to answer is whether there is a phase transition (first order). This is unresolved currently. Another question we might ask is whether there exists multiple ground state pairs (i.e. we can find 2 different configurations that have minimal energy that are not merely negative images of each other) - again this is unresolved. In infinite dimensions we can see this would be true, in 1d we know this cannot be true - for other dimensions it is unclear (although it is believed it is probably not true for d=2).
In addition to this there are many other unanswered questions that are the subject of current research into Ising type models.
Conclusion
In this blog post we have investigated some of the properties of square-lattice Ising models in various dimensions. In particular we have seen that there is no phase transition in 1d, a second order phase tranisition in 2d, a phase transition in infinite dimension and currently other dimensions are unresolved. We can see that the short-range interactions cause a lot of headaches when trying to analyse these systems. In the next blog post in this series we will begin to look at ways of simulating Ising models (and other spin glass models generally).
This is the third blog post in a series - you can find the next blog post here