The Ising model

Ferromagnetism arises when a collection of atomic spins align such that their associated magnetic moments all point in the same direction, yielding a net magnetic moment which is macroscopic in size. The simplest theoretical description of ferromagnetism is called the Ising model. This model was invented by Wilhelm Lenz in 1920: it is named after Ernst Ising, a student of Lenz who chose the model as the subject of his doctoral dissertation in 1925.

Consider $N$ atoms in the presence of a $z$-directed magnetic field of strength $H$. Suppose that all atoms are identical spin-$1/2$ systems. It follows that either $s_i=+1$ (spin up) or $s_i=-1$ (spin down), where $s_i$ is (twice) the $z$-component of the $i$th atomic spin. The total energy of the system is written:

$$E = - J\sum_{<i j>}s_i\,s_j -\mu\,H\sum_{i=1,N}s_i.$$ (351)

Here, $<i j>$ refers to a sum over nearest neighbour pairs of atoms. Furthermore, $J$ is called the exchange energy, whereas $\mu$ is the atomic magnetic moment. Equation (351) is the essence of the Ising model.

The physics of the Ising model is as follows. The first term on the right-hand side of Eq. (351) shows that the overall energy is lowered when neighbouring atomic spins are aligned. This effect is mostly due to the Pauli exclusion principle. Electrons cannot occupy the same quantum state, so two electrons on neighbouring atoms which have parallel spins (i.e., occupy the same orbital state) cannot come close together in space. No such restriction applies if the electrons have anti-parallel spins. Different spatial separations imply different electrostatic interaction energies, and the exchange energy, $J$, measures this difference. Note that since the exchange energy is electrostatic in origin, it can be quite large: i.e., $J\sim 1$eV. This is far larger than the energy associated with the direct magnetic interaction between neighbouring atomic spins, which is only about $10^{-4}$eV. However, the exchange effect is very short-range; hence, the restriction to nearest neighbour interaction is quite realistic.

Our first attempt to analyze the Ising model will employ a simplification known as the mean field approximation. The energy of the $i$th atom is written

$$e_i = -\frac{J}{2}\sum_{k=1,z}s_k\,s_i -\mu\,H\,s_i,$$ (352)

where the sum is over the $z$ nearest neighbours of atom $i$. The factor $1/2$ is needed to ensure that when we sum to obtain the total energy,

$$E =\sum_{i=1,N} e_i,$$ (353)

we do not count each pair of neighbouring atoms twice.

We can write

$$e_i = -\mu\,H_{\rm eff}\,s_i,$$ (354)

where

$$H_{\rm eff} = H + \frac{J}{2\,\mu}\sum_{k=1,z}s_k.$$ (355)

Here, $H_{\rm eff}$ is the effective magnetic field, which is made up of two components: the external field, $H$, and the internal field generated by neighbouring atoms.

Consider a single atom in a magnetic field $H_m$. Suppose that the atom is in thermal equilibrium with a heat bath of temperature $T$. According to the well-known Boltzmann distribution, the mean spin of the atom is

$$\bar{s } = \frac{{\rm e}^{+\beta\,\mu\,H_m} - {\rm e}^{-\bet... ...m}} {{\rm e}^{+\beta\,\mu\,H_m} + {\rm e}^{-\beta\,\mu\,H_m}},$$ (356)

where $\beta = 1/kT$, and $k$ is the Boltzmann constant. The above expression follows because the energy of the ``spin up'' state ($s=+1$) is $-\mu\,H_m$, whereas the energy of the ``spin down'' state ($s=-1$) is $+\mu\,H_m$. Hence,

$$\bar{s} = \tanh(\beta\,\mu\,H_m).$$ (357)

Let us assume that all atoms have identical spins: i.e., $s_i=\bar{s}$. This assumption is known as the ``mean field approximation''. We can write

$$H_{\rm eff} = H + \frac{z\,J\,\bar{s}}{2\,\mu}.$$ (358)

Finally, we can combine Eqs. (357) and (358) (identifying $H_m$ and $H_{\rm eff}$) to obtain

$$\bar{s} = \tanh\left\{\beta\,\mu\,H + \beta\,z\,J\,\bar{s}/2\right\}.$$ (359)

Note that the heat bath in which a given atom is immersed is simply the rest of the atoms. Hence, $T$ is the temperature of the atomic array. It is helpful to define the critical temperature,

$$T_c = \frac{z\,J}{2\,k},$$ (360)

and the critical magnetic field,

$$H_c = \frac{k\,T_c}{\mu} = \frac{z\,J}{2\,\mu}.$$ (361)

Equation (359) reduces to

$$\bar{s} = \tanh\left\{\frac{T_c}{T}\left(\frac{H}{H_c} + \bar{s}\right)\right\}.$$ (362)

The above equation cannot be solved analytically. However, it is fairly easily to solve numerically using the following iteration scheme:

$$\bar{s}_{i+1} = \tanh\left\{\frac{T_c}{T}\left(\frac{H}{H_c} + \bar{s}_i\right)\right\}.$$ (363)

The above formula is iterated until $\bar{s}_{i+1}\rightarrow \bar{s}_i$.

It is helpful to define the net magnetization,

$$M=\mu\sum_{i=1,N} s_i= \mu\,N\,\bar{s},$$ (364)

the net energy,

$$E=\sum_{i=1,N} e_i= -N\,k\,T_c\left(\frac{H}{H_c}+\bar{s}\right)\bar{s},$$ (365)

and the heat capacity,

$$C=\frac{dE}{dT}.$$ (366)

Figure 101: The net magnetization, $M$, of a collection of $N$ ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Calculation performed using the mean field approximation.

Figures 102, 101, and 103 show the net magnetization, net energy, and heat capacity calculated from the iteration formula (363) in the absence of an external magnetic field (i.e., with $H=0$). It can be seen that below the critical (or ``Curie'') temperature, $T_c$, there is spontaneous magnetization: i.e., the exchange effect is sufficiently large to cause neighbouring atomic spins to spontaneously align. On the other hand, thermal fluctuations completely eliminate any alignment above the critical temperature. Moreover, at the critical temperature there is a discontinuity in the first derivative of the energy, $E$, with respect to the temperature, $T$. This discontinuity generates a downward jump in the heat capacity, $C$, at $T=T_c$. The sudden loss of spontaneous magnetization as the temperature exceeds the critical temperature is a type of phase transition.

Figure 102: The net energy, $E$, of a collection of $N$ ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Calculation performed using the mean field approximation.

Now, according to the conventional classification of phase transitions, a transition is first-order if the energy is discontinuous with respect to the order parameter (i.e., in this case, the temperature), and second-order if the energy is continuous, but its first derivative with respect to the order parameter is discontinuous, etc. We conclude that the loss of spontaneous magnetization in a ferromagnetic material as the temperature exceeds the critical temperature is a second-order phase transition.

Figure 103: The heat capacity, $C$, of a collection of $N$ ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Calculation performed using the mean field approximation.

In order to see an example of a first-order phase transition, let us examine the behaviour of the magnetization, $M$, as the external field, $H$, is varied at constant temperature, $T$.

Figure 104: The net magnetization, $M$, of a collection of $N$ ferromagnetic atoms as a function of the external magnetic field, $H$, at constant temperature, $T<T_c$. Calculation performed using the mean field approximation.

Figure 105: The net energy, $E$, of a collection of $N$ ferromagnetic atoms as a function of the external magnetic field, $H$, at constant temperature, $T<T_c$. Calculation performed using the mean field approximation.

Figures 104 and 105 show the magnetization, $M$, and energy, $E$, versus external field-strength, $H$, calculated from the iteration formula (363) at some constant temperature, $T$, which is less than the critical temperature, $T_c$. It can be seen that $E$ is discontinuous, indicating the presence of a first-order phase transition. Moreover, the system exhibits hysteresis--meta-stable states exist within a certain range of $H$ values, and the magnetization of the system at fixed $T$ and $H$ (within the aforementioned range) depends on its past history: i.e., on whether $H$ was increasing or decreasing when it entered the meta-stable range.

Figure 106: The net magnetization, $M$, of a collection of $N$ ferromagnetic atoms as a function of the external magnetic field, $H$, at constant temperature, $T=T_c$. Calculation performed using the mean field approximation.

Figure 107: The net energy, $E$, of a collection of $N$ ferromagnetic atoms as a function of the external magnetic field, $H$, at constant temperature, $T=T_c$. Calculation performed using the mean field approximation.

Figures 106 and 107 show the magnetization, $M$, and energy, $E$, versus external field-strength, $H$, calculated from the iteration formula (363) at a constant temperature, $T$, which is equal to the critical temperature, $T_c$. It can be seen that $E$ is now continuous, and there are no meta-stable states. We conclude that first-order phase transitions and hysteresis only occur, as the external field-strength is varied, when the temperature lies below the critical temperature: i.e., when the ferromagnetic material in question is capable of spontaneous magnetization.

The above calculations, which are based on the mean field approximation, correctly predict the existence of first- and second-order phase transitions when $H\neq 0$ and $H=0$, respectively. However, these calculations get some of the details of the second-order phase transition wrong. In order to do a better job, we must abandon the mean field approximation and adopt a Monte-Carlo approach.

Figure 108: A two-dimensional array of atoms.

Let us consider a two-dimensional square array of atoms. Let $L$ be the size of the array, and $N=L^2$ the number of atoms in the array, as shown in Fig. 108. The Monte-Carlo approach to the Ising model, which completely avoids the use of the mean field approximation, is based on the following algorithm:

• Step through each atom in the array in turn:
  • For a given atom, evaluate the change in energy of the system, ${\mit\Delta}E$, when the atomic spin is flipped.
  • If ${\mit\Delta}E<0$ then flip the spin.
  • If ${\mit\Delta}E>0$ then flip the spin with probability $P=\exp(-\beta\,{\mit\Delta}E)$.
• Repeat the process many times until thermal equilibrium is achieved.

The purpose of the algorithm is to shuffle through all possible states of the system, and to ensure that the system occupies a given state with the Boltzmann probability: i.e., with a probability proportional to $\exp(-\beta\,E)$, where $E$ is the energy of the state.

In order to demonstrate that the above algorithm is correct, let us consider flipping the spin of the $i$th atom. Suppose that this operation causes the system to make a transition from state $a$ (energy, $E_a$) to state $b$ (energy, $E_b$). Suppose, further, that $E_a < E_b$. According to the above algorithm, the probability of a transition from state $a$ to state $b$ is

$$P_{a\rightarrow b} = \exp[-\beta\,(E_b-E_a)],$$ (367)

whereas the probability of a transition from state $b$ to state $a$ is

$$P_{b\rightarrow a} = 1.$$ (368)

In thermal equilibrium, the well-known principal of detailed balance implies that

$$P_a\,P_{a\rightarrow b} = P_b\,P_{b\rightarrow a},$$ (369)

where $P_a$ is the probability that the system occupies state $a$, and $P_b$ is the probability that the system occupies state $b$. Equation (369) simply states that in thermal equilibrium the rate at which the system makes transitions from state $a$ to state $b$ is equal to the rate at which the system makes reverse transitions. The previous equation can be rearranged to give

$$\frac{P_b}{P_a} = \exp[-\beta\,(E_b-E_a)],$$ (370)

which is consistent with the Boltzmann distribution.

Now, each atom in our array has four nearest neighbours, except for atoms on the edge of the array, which have less than four neighbours. We can eliminate this annoying special behaviour by adopting periodic boundary conditions: i.e., by identifying opposite edges of the array. Indeed, we can think of the array as existing on the surface of a torus.

It is helpful to define

$$T_0 = \frac{J}{k}.$$ (371)

Now, according to mean field theory,

$$T_c = \frac{z\,J}{2\,k}= 2\,T_0.$$ (372)

The evaluation of

$$C = \lim_{{\mit\Delta}T\rightarrow 0}\frac{{\mit\Delta}E}{{\mit\Delta} T}$$ (373)

via the direct method is difficult due to statistical noise in the energy, $E$. Instead, we can make use of a standard result in equilibrium statistical thermodynamics:

$$C = \frac{\sigma_{E}^{\,2}}{k\,T^2},$$ (374)

where $\sigma_E$ is the standard deviation of fluctuations in $E$. Fortunately, it is fairly easy to evaluate $\sigma_E$: we can simply employ the standard deviation in $E$ from step to step in our Monte-Carlo iteration scheme.

Figure 109: The net magnetization, $M$, of a $5\times 5$ array of ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Monte-Carlo simulation.

Figures 109-116 show magnetization and heat capacity versus temperature curves for $L=5$, 10, 20, and 40 in the absence of an external magnetic field. In all cases, the Monte-Carlo simulation is iterated 5000 times, and the first 1000 iterations are discarded when evaluating $\sigma_E$ (in order to allow the system to attain thermal equilibrium). The two-dimensional array of atoms is initialized in a fully aligned state for each different value of the temperature. Since there is no external magnetic field, it is irrelevant whether the magnetization, M, is positive or negative. Hence, $M$ is replaced by $\vert M\vert$ in all plots.

Figure 110: The heat capacity, $C$, of a $5\times 5$ array of ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Monte-Carlo simulation. The solid curve shows the heat capacity calculated from Eq. (373), whereas the dotted curve shows the heat capacity calculated from Eq. (374).

Note that the $M$ versus $T$ curves generated by the Monte-Carlo simulations look very much like those predicted by the mean field model. The resemblance increases as the size, $L$, of the atomic array increases. The major difference is the presence of a magnetization ``tail'' for $T>T_c$ in the Monte-Carlo simulations: i.e., in the Monte-Carlo simulations the spontaneous magnetization does not collapse to zero once the critical temperature is exceeded--there is a small lingering magnetization for $T>T_c$. The $C$ versus $T$ curves show the heat capacity calculated directly (i.e., $C={\mit\Delta E}/{\mit\Delta T}$), and via the identity $C=\sigma_E^{\,2}/k\,T^2$. The latter method of calculation is clearly far superior, since it generates significantly less statistical noise. Note that the heat capacity peaks at the critical temperature: i.e., unlike the mean field model, $C$ is not zero for $T>T_c$. This effect is due to the residual magnetization present when $T>T_c$.

Figure 111: The net magnetization, $M$, of a $10\times 10$ array of ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Monte-Carlo simulation.

Our best estimate for $T_c$ is obtained from the location of the peak in the $C$ versus $T$ curve in Fig. 116. We obtain $T_c=2.27\,T_0$. Recall that the mean field model yields $T_c = 2\,T_0$. The exact answer for a two-dimensional array of ferromagnetic atoms is

$$T_c = \frac{2\,T_0}{\ln(1+\sqrt{2})} = 2.27\,T_0,$$ (375)

which is consistent with our Monte-Carlo calculations. The above analytic result was first obtained by Onsager in 1944.³⁹ Incidentally, Onsager's analytic solution of the 2-D Ising model is one of the most complicated and involved calculations in all of theoretical physics. Needless to say, no one has ever been able to find an analytic solution of the Ising model in more than two dimensions.

Figure 112: The heat capacity, $C$, of a $10\times 10$ array of ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Monte-Carlo simulation. The solid curve shows the heat capacity calculated from Eq. (373), whereas the dotted curve shows the heat capacity calculated from Eq. (374).

Note, from Figs. 110, 112, 114, and 116, that the height of the peak in the heat capacity curve at $T=T_c$ increases with increasing array size, $L$. Indeed, a close examination of these figures yields $C_{\rm max}/N\,k = 0.95$ for $L=5$, $C_{\rm max}/N\,k = 1.34$ for $L=10$, $C_{\rm max}/N\,k = 1.77$ for $L=20$, and $C_{\rm max}/N\,k = 2.16$ for $L=40$. Figure 117 shows $C_{\rm max}/N\,k$ plotted against $\ln L$ for $L=5$, 10, 20, and 40. It can be seen that the points lie on a very convincing straight-line, which strongly suggests that

$$\frac{C_{\rm max}}{k\,N} \propto \ln L.$$ (376)

Figure 113: The net magnetization, $M$, of a $20\times 20$ array of ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Monte-Carlo simulation.

Of course, for physical systems, $L\sim \sqrt{N_A} \sim 10^{12}$, where $N_A$ is Avogadro's number. Hence, $C$ is effectively singular at the critical temperature (since $\ln N_A \gg 1$), as sketched in Fig. 118. This observation leads us to revise our definition of a second-order phase transition. It turns out that actual discontinuities in the heat capacity almost never occur. Instead, second-order phase transitions are characterized by a local quasi-singularity in the heat capacity.

Figure 114: The heat capacity, $C$, of a $20\times 20$ array of ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Monte-Carlo simulation. The solid curve shows the heat capacity calculated from Eq. (373), whereas the dotted curve shows the heat capacity calculated from Eq. (374).

Recall, from Eq. (374), that the typical amplitude of energy fluctuations is proportional to the square-root of the heat capacity (i.e., $\sigma_E\propto \sqrt{C}$). It follows that the amplitude of energy fluctuations becomes extremely large in the vicinity of a second-order phase transition.

Figure 115: The net magnetization, $M$, of a $40\times 40$ array of ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Monte-Carlo simulation.

Now, the main difference between our mean field and Monte-Carlo calculations is the existence of residual magnetization for $T>T_c$ in the latter case. Figures 119-123 show the magnetization pattern of a $40\times 40$ array of ferromagnetic atoms, in thermal equilibrium and in the absence of an external magnetic field, calculated at various temperatures. It can be seen that for $T=20\,T_0$ the pattern is essentially random. However, for $T=5\,T_0$, small clumps appear in the pattern. For $T=3\,T_0$, the clumps are somewhat bigger. For $T=2.32\,T_0$, which is just above the critical temperature, the clumps are global in extent. Finally, for $T=1.8\,T_0$, which is a little below the critical temperature, there is almost complete alignment of the atomic spins.

Figure 116: The heat capacity, $C$, of a $40\times 40$ array of ferromagnetic atoms as a function of the temperature, $T$, in the absence of an external magnetic field. Monte-Carlo simulation. The solid curve shows the heat capacity calculated from Eq. (373), whereas the dotted curve shows the heat capacity calculated from Eq. (374).

The problem with the mean field model is that it assumes that all atoms are situated in identical environments. Hence, if the exchange effect is not sufficiently large to cause global alignment of the atomic spins then there is no alignment at all. What actually happens when the temperature exceeds the critical temperature is that global alignment disappears, but local alignment (i.e., clumping) remains. Clumps are only eliminated by thermal fluctuations once the temperature is significantly greater than the critical temperature. Atoms in the middle of the clumps are situated in a different environment than atoms on the clump boundaries. Hence, clumps cannot occur in the mean field model.

Figure 117: The peak value of the heat capacity (normalized by $N\,k$) versus the logarithm of the array size for a two-dimensional array of ferromagnetic atoms in the absence of an external magnetic field. Monte-Carlo simulation.

Figure 118: A sketch of the expected variation of the heat capacity versus the temperature for a physical two-dimensional ferromagnetic system.

Figure 119: Magnetization pattern of a $40\times 40$ array of ferromagnetic atoms in thermal equilibrium and in the absence of an external magnetic field. Monte-Carlo calculation with $T=20\,T_0$. Black/white squares indicate atoms magnetized in plus/minus $z$-direction, respectively.

Figure 120: Magnetization pattern of a $40\times 40$ array of ferromagnetic atoms in thermal equilibrium and in the absence of an external magnetic field. Monte-Carlo calculation with $T=5\,T_0$. Black/white squares indicate atoms magnetized in plus/minus $z$-direction, respectively.

Figure 121: Magnetization pattern of a $40\times 40$ array of ferromagnetic atoms in thermal equilibrium and in the absence of an external magnetic field. Monte-Carlo calculation with $T=3\,T_0$. Black/white squares indicate atoms magnetized in plus/minus $z$-direction, respectively.

Figure 122: Magnetization pattern of a $40\times 40$ array of ferromagnetic atoms in thermal equilibrium and in the absence of an external magnetic field. Monte-Carlo calculation with $T=2.32\,T_0$. Black/white squares indicate atoms magnetized in plus/minus $z$-direction, respectively.

Figure 123: Magnetization pattern of a $40\times 40$ array of ferromagnetic atoms in thermal equilibrium and in the absence of an external magnetic field. Monte-Carlo calculation with $T=1.8\,T_0$. Black/white squares indicate atoms magnetized in plus/minus $z$-direction, respectively.