Project 4

First and second order phase transitions in the Ising model of ferromagnetism

Introduction

Write a program which simulates via Metropolis Monte Carlo algorithm the 2D Ising model on a $128\times 128$ square lattice with periodic boundary conditions. The model is assumed to describe a system placed in an applied external magnetic field $H$. The energy of a particular state of this system is then given by:

$E=-J\sum _{\langle i,j\rangle }S_{i}S_{j}-H\sum _{i}S_{i}$

where $\langle i,j\rangle$ means that the sum runs over all pairs of nearest neighbor spins (attached to the lattice sites), and $J$ is the exchange coupling. Positive $J>0$ interaction between the spins leads to ferromagnetism, will negative $J>0$ gives rise to antiferromagnetism in the Ising model. Choose units so that $|J|=1$ and the Boltzmann constant is $k_{B}=1$. The free parameters of the model are then just the temperature $T$ and the applied field $H$.

Part I for both PHYS460 and PHYS660 students: First order phase transition in the ferromagnetic Ising model

Consider the ferromagnetic Ising model with $J=1$. Start from a case in which temperature is below the critical one $T_{c}=2.269$, such as $T=1.0$, and the field is large and pointing down, e.g., $H=-5$. The selected initial state of the Ising model should be the one in which all spins $S_{i}=-1$ are parallel to H, so that magnetization per spin is $M=-1$.

Sweep the magnetic field up to a large positive value $H=5$ and then back down to your starting value $H=-5$. Do this in steps of $\delta H=0.5$. After each change in the magnetic field $H$, allow sufficient Monte Carlo trials for the system to reach the equilibrium and thermalize.

For each value of $H$, the mean magnetization per spin

$\langle M\rangle ={\frac {1}{N_{MC}}}\sum _{\alpha =1}^{N_{MC}}M_{\alpha }$

The sum in this expression run are over $N_{MC}$ microstates generated by Monte Carlo simulation.

(a) Plot $M$ vs. $H$. You should see a discontinuous change in as H is increased, and a second discontinuous change in $M$ as $H$ is decreased back to its initial value. These discontinuous changes indicate a first order phase transition. Also note that the values of $H$ at which discontinuities occur are not identical. Such behavior is called hysteresis. Perform at least one thousand Monte Carlo time steps, where one time step is one complete pass through the lattice, at each value of the magnetic field.

(b) Repeat the calculations for a temperature $T=4.0$ above the critical point $T_{c}$. Is there a first order phase transition? Is there hysteresis? Provide a physical explanation of any difference in behavior for the two temperatures.

Part II for both PHYS460 and PHYS660 students: Second order phase transition in the ferromagnetic Ising model

(a) Plot magnetization per spin and energy per spin for each sample obtained after one Monte Carlo sweep through the lattice as a function of Monte Carlo time step at three temperatures $T<T_{c}$, $T=T_{c}$, and $T>T_{c}$ in the absence of magnetic field $H=0$. These plots should reveal how many Monte Carlo steps is necessary to thermalize (i.e., bring it into the thermodynamic equilibrium with external heat bath which is assumed to set the temperature $T$) the Ising ferromagnet. Thermodynamic average values should be computed only for the samples generated after the thermalization has been achieved.

(b) Investigate the behavior of the specific heat

$C_{V}={\frac {1}{k_{B}T^{2}}}(\langle E^{2}\rangle -\langle E\rangle ^{2})$

and magnetic susceptibility

$\chi ={\frac {1}{k_{B}T}}(\langle M^{2}\rangle -\langle M\rangle ^{2})$

as the function of increasing temperature in the absence of external magnetic field $H=0$. In this case we expect a second order or continuous paramagnet-ferromanget phase transition at the critical temperature $T_{c}$ around which these two quantities diverge in the thermodynamic limit. Since your calculations is for finite size system, repeat computation of these two quantities while increasing the lattice size (e.g., $16\times 16$, $32\times 32$, $64\times 64$, ...) in order to observe how the peak in $C_{V}$ and $\chi$ around $T_{c}$ increases while its center shifts because toward $T_{c}$. The two quantities should be divided by the lattice size to obtain their values "per spin".

The connection between these two response functions and fluctuations

$(\Delta X)^{2}=\langle X^{2}\rangle -\langle X\rangle ^{2}$

where

$\langle X^{2}\rangle ={\frac {1}{N_{MC}}}\sum _{\alpha =1}^{N_{MC}}X_{\alpha }^{2}$

$\langle X\rangle ^{2}=\left({\frac {1}{N_{MC}}}\sum _{\alpha =1}^{N_{MC}}X_{\alpha }\right)^{2}$

utilizes the so-called fluctuation-dissipation theorem of statistical physics.

Part III for PHYS660 students only: Correlation length in the ferromagnetic Ising model

The correlation lenght can be obtained from the r-dependence of the spin-spin correlation function

$G(r)=\langle S_{i}S_{j}\rangle -\langle S_{i}\rangle \langle S_{j}\rangle$

where $\langle S_{i}\rangle \langle S_{j}\rangle =M^{2}$ since the system is translationally invariant. Here $r=|i-j|$ is the distance between the sites $i$ and $j$. The average here $\langle \ldots \rangle$ is defined over all sites for a given configuration and over many configurations. Because the spins are not correlated for large $r$, $G(r)\rightarrow 0$ in this limit. Assume that $G(r)\sim \exp(-r/\xi )$ for sufficiently large $r$ and estimate as a function of temperature $T$ in zero external field $H$.

How does your estimate for $\xi$ compare with the size of the domains of spins pointing in the same direction? We expect that $\xi$ should be of the order of a lattice spacing at $T\gg T_{c}$ and approaching the system size (or diverging for infinite systems in the thermodynamic limit) as $T\rightarrow T_{c}$.

Part IV for PHYS660 students only: Second order phase transition in the antiferromagnetic Ising model

Compute temperature dependence of $C_{V}$, $\chi$, and $M$ in the absence of magnetic field $H=0$ for the antiferromagnetic Ising model with $J=-1$. The initial configuration can be chosen as all spins up $S_{i}=+1$ on a $128\times 128$ lattice. What configuration of spins corresponds to the thermodynamic state below the Neel temperature $T_{N}=2.269$? At the so-called Neel temperature this type of materials exhibit continous antiferromagnet-paramagnet phase transition.

References

• I. Vattulainen, T. Ala-Nissila, and K. Kankaala, Physical tests for random numbers in simulations, Phys. Rev. Lett. 73, 2513 (1994). [PDF]
• D. P. Landau, S.-H. Tsai, and M. Exler, A new approach to Monte Carlo simulations in statistical physics: Wang-Landau sampling, Am. J. Phys. 72, 1294 (2004). [PDF]