Ising Model

Statistical Mechanics

The Ising model is a mathematical model of ferromagnetism in statistical mechanics. It consists of discrete variables (spins) that can be in one of two states (+1 or -1) arranged on a lattice, with interactions between neighboring spins.

Hamiltonian

For a 2D square lattice, the Ising model Hamiltonian is:

H = - J ⟨ i, j ⟩ \sum s_{i} s_{j} - h i \sum s_{i}

where:

$J$ is the coupling constant (ferromagnetic if $J > 0$ )
$s_{i} \in {- 1, + 1}$ are the spin variables
$⟨ i, j ⟩$ denotes nearest neighbor pairs
$h$ is the external magnetic field

Metropolis Algorithm

The Metropolis algorithm is used to sample spin configurations according to the Boltzmann distribution. At each step:

Select a random spin
Calculate the energy change $Δ E$ if the spin is flipped
Accept the flip with probability $min (1, e^{- β Δ E})$ where $β = 1/ k_{B} T$

The energy change for flipping spin $s_{i}$ is:

Δ E = 2 s_{i} J j \in neighbors \sum s_{j} + h

Phase Transition

The 2D Ising model exhibits a phase transition at the critical temperature $T_{c} \approx 2.269 J / k_{B}$ . Below $T_{c}$ , the system is ferromagnetic (spontaneous magnetization), and above $T_{c}$ , it is paramagnetic (no net magnetization).

Binder Cumulant

The Binder cumulant is a useful quantity for identifying phase transitions:

U_{L} = 1 - \frac{⟨ m ^{4} ⟩ _{L}}{3 ⟨ m ^{2} ⟩ _{L}^{2}}

where $m$ is the magnetization per spin and $L$ is the system size. At the critical temperature, the Binder cumulant becomes size-independent, making it useful for finite-size scaling analysis.

Python Implementation

The following code implements the 2D Ising model with Metropolis Monte Carlo sampling:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

def set_publication_style():
    """Set publication-quality matplotlib style."""
    mpl.rcParams.update({
        'font.family': 'serif',
        'font.size': 12,
        'axes.labelsize': 14,
        'axes.titlesize': 16,
        'axes.linewidth': 1.2,
        'axes.labelpad': 8,
        'axes.titlepad': 10,
        'xtick.labelsize': 12,
        'ytick.labelsize': 12,
        'xtick.direction': 'in',
        'ytick.direction': 'in',
        'xtick.top': True,
        'ytick.right': True,
        'xtick.major.size': 6,
        'ytick.major.size': 6,
        'xtick.major.width': 1.2,
        'ytick.major.width': 1.2,
        'legend.fontsize': 12,
        'legend.frameon': False,
        'lines.linewidth': 2,
        'lines.markersize': 6,
        'figure.dpi': 100,
        'savefig.dpi': 300,
        'savefig.bbox': 'tight'
    })

set_publication_style()

def initialize_spins(N):
    """Initialize a random spin configuration."""
    return np.random.choice([-1, 1], size=(N, N))

def calculate_magnetization(spins):
    """Calculate total magnetization."""
    return np.sum(spins)

def calculate_energy(spins, J=1.0, h=0.0):
    """Calculate the energy of a spin configuration."""
    N = spins.shape[0]
    energy = 0.0
    
    # Nearest neighbor interactions (periodic boundary conditions)
    for i in range(N):
        for j in range(N):
            energy -= J * spins[i, j] * (
                spins[(i+1) % N, j] +  # Right neighbor
                spins[(i-1) % N, j] +  # Left neighbor
                spins[i, (j+1) % N] +  # Top neighbor
                spins[i, (j-1) % N]    # Bottom neighbor
            )
    
    # External field contribution
    energy -= h * np.sum(spins)
    
    return energy / 2  # Each pair counted twice

def metropolis_step(spins, beta, J=1.0, h=0.0):
    """Perform one Metropolis Monte Carlo step."""
    N = spins.shape[0]
    
    for _ in range(N * N):  # One sweep through all spins
        i, j = np.random.randint(0, N, size=2)
        
        # Calculate energy change if spin flips
        spin_flip_energy = 2 * spins[i, j] * (
            J * (spins[(i+1) % N, j] + spins[(i-1) % N, j] +
                 spins[i, (j+1) % N] + spins[i, (j-1) % N]) + h
        )
        
        # Metropolis acceptance criterion
        if spin_flip_energy <= 0 or np.random.rand() < np.exp(-spin_flip_energy * beta):
            spins[i, j] *= -1
    
    return spins

def simulate_ising(N, T, num_steps, thermalization_steps, J=1.0, h=0.0):
    """Simulate the Ising model and collect statistics."""
    spins = initialize_spins(N)
    beta = 1.0 / T
    
    magnetizations = []
    magnetizations_squared = []
    magnetizations_fourth = []
    
    for step in range(num_steps):
        spins = metropolis_step(spins, beta, J, h)
        
        if step >= thermalization_steps:
            m = calculate_magnetization(spins) / (N * N)  # Magnetization per spin
            magnetizations.append(m)
            magnetizations_squared.append(m**2)
            magnetizations_fourth.append(m**4)
    
    return (np.array(magnetizations),
            np.array(magnetizations_squared),
            np.array(magnetizations_fourth))

# Parameters
N = 20  # Lattice size (N x N)
temperatures = np.linspace(1.0, 3.5, 30)  # Temperature range
num_steps = 20000  # Total Monte Carlo steps
thermalization_steps = 5000  # Steps to discard for thermalization
J = 1.0  # Coupling constant
h = 0.0  # External field

# Arrays to store results
mean_magnetization = []
susceptibility = []
binder_cumulant = []

# Simulate at different temperatures
for T in temperatures:
    m, m2, m4 = simulate_ising(N, T, num_steps, thermalization_steps, J, h)
    
    mean_m = np.mean(m)
    mean_m2 = np.mean(m2)
    mean_m4 = np.mean(m4)
    
    # Susceptibility: chi = (1/T) * (<m^2> - <m>^2)
    chi = (mean_m2 - mean_m**2) / T
    
    # Binder cumulant
    U = 1 - mean_m4 / (3 * mean_m2**2) if mean_m2 > 0 else 0
    
    mean_magnetization.append(np.abs(mean_m))  # Use absolute value for magnetization
    susceptibility.append(chi)
    binder_cumulant.append(U)

# Plotting
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Magnetization vs Temperature
axes[0].plot(temperatures, mean_magnetization, 'o-', linewidth=2, markersize=4)
axes[0].axvline(x=2.269, color='r', linestyle='--', label='T_c ≈ 2.269')
axes[0].set_xlabel('Temperature (T)')
axes[0].set_ylabel('|Magnetization| per Spin')
axes[0].set_title('Magnetization vs Temperature')
axes[0].legend()
axes[0].grid(True)

# Susceptibility vs Temperature
axes[1].plot(temperatures, susceptibility, 's-', linewidth=2, markersize=4, color='green')
axes[1].axvline(x=2.269, color='r', linestyle='--', label='T_c ≈ 2.269')
axes[1].set_xlabel('Temperature (T)')
axes[1].set_ylabel('Susceptibility (χ)')
axes[1].set_title('Susceptibility vs Temperature')
axes[1].legend()
axes[1].grid(True)

# Binder Cumulant vs Temperature
axes[2].plot(temperatures, binder_cumulant, '^-', linewidth=2, markersize=4, color='purple')
axes[2].axvline(x=2.269, color='r', linestyle='--', label='T_c ≈ 2.269')
axes[2].set_xlabel('Temperature (T)')
axes[2].set_ylabel('Binder Cumulant (U)')
axes[2].set_title('Binder Cumulant vs Temperature')
axes[2].legend()
axes[2].grid(True)

plt.tight_layout()
plt.savefig('figures/ising_model_phase_transition.png', dpi=300, bbox_inches='tight')
plt.show()

print(f"Critical temperature (exact): T_c = 2.269 J/k_B")
print(f"Peak susceptibility at T ≈ {temperatures[np.argmax(susceptibility)]:.3f}")

Visualization

The following plot shows the phase transition behavior of the Ising model:

Figure pending — running the script above produces figures/ising_model_phase_transition.png.

Key Features

Metropolis Monte Carlo algorithm for sampling
Periodic boundary conditions
Thermalization period to reach equilibrium
Calculation of magnetization, susceptibility, and Binder cumulant
Phase transition detection

Physical Significance

The Ising model is one of the simplest models that exhibits a phase transition. It has applications in:

Magnetic materials
Lattice gas models
Binary alloys
Neural networks (Hopfield model)
Image processing

Critical Exponents

Near the critical temperature, various quantities exhibit power-law behavior:

Magnetization: $M \sim ∣ T - T_{c} ∣^{β}$ with $β = 1/8$ (2D)
Susceptibility: $χ \sim ∣ T - T_{c} ∣^{- γ}$ with $γ = 7/4$ (2D)
Correlation length: $ξ \sim ∣ T - T_{c} ∣^{- ν}$ with $ν = 1$ (2D)