Introduction

We simulate the 3D Ising model on an N x N x N cubic lattice using the Metropolis Monte Carlo algorithm. We compute the magnetization per spin as a function of temperature for different lattice sizes as the main observable to monitor, which will make manifest the phase transition which is the characteristic feature of this model.

# Imports

import numpy as np
import math
import matplotlib.pyplot as plt
from numba import jit, njit
from joblib import Parallel, delayed
from tqdm.notebook import tqdm

SEED = 33
np.random.seed(SEED)

We define some functions to compute quantities of interest, namely:

$\cdot$ the energy of the configuration, given by the Hamiltonian: $$H=-J\sum_{\langle i, j\rangle}S_iS_j$$ where we sum on neighboring 'neurons';

$\cdot$ the difference in energy $\Delta E$ when flipping a single spin: i.e. $S_i \longrightarrow -S_i$. Notice that, after a single flip spin, we return the difference in energy $$\Delta E=2JS_i\sum_{\langle i, j\rangle}S_j$$ in fact, before: $E_{\text{init.}}=-JS_i\sum_{\langle i, j\rangle}S_j$ and after: $E_{\text{fin.}}=-J(-S_i)\sum_{\langle i, j\rangle}S_j$, hence the expression;

$\cdot$ the average magnetization of each configuration, namely, the mean of the total spin sum;

$\cdot$ the autocorrelation function, i.e., the weighted covariance between the atoms.

# JIT-compiled energy difference and update

@jit(nopython=True)
def compute_energy(lattice, J):
    E = 0.0
    N = lattice.shape[0]
    for x in range(N):
        for y in range(N):
            for z in range(N):
                S = lattice[x,y,z]
                E -= J * S * (
                    lattice[(x+1)%N,y,z] + lattice[x,(y+1)%N,z] + lattice[x,y,(z+1)%N]
                )
    return E

@jit(nopython=True)
def delta_energy(lattice, x, y, z, J):
    N = lattice.shape[0]
    S = lattice[x, y, z]
    neighbors = (
        lattice[(x+1)%N, y, z] + lattice[(x-1)%N, y, z] +
        lattice[x, (y+1)%N, z] + lattice[x, (y-1)%N, z] +
        lattice[x, y, (z+1)%N] + lattice[x, y, (z-1)%N]
    )
    return 2 * J * S * neighbors

@jit(nopython=True)
def compute_magnetization(lattice):
    return np.abs(lattice.sum()) / lattice.size

def autocorrelation(data):
    n = len(data)
    data = np.array(data) - np.mean(data)
    var = np.var(data)
    acf = np.correlate(data, data, mode='full')[n-1:] / (var * np.arange(n, 0, -1))
    return acf

The following function does a single sweep of the Metropolis algorithm. At each iteration it picks a spin uniformly at random and looks at the change in energy. The spin is then flipped according to the Metropolis rule:

$\cdot$ if the energy is decreased, i.e. $\Delta E<0$, we always accept the proposed move;

$\cdot$, we accept the move with probability $\exp(-\Delta E/T)$.

@njit
def metropolis_sweep(lattice, beta, J, xs, ys, zs, us):
    N = lattice.shape[0]
    for i in range(N**3):
        x = xs[i]
        y = ys[i]
        z = zs[i]
        # compute ΔE via jit‐ted delta_energy
        dE = delta_energy(lattice, x, y, z, J)
        if dE <= 0 or us[i] < math.exp(-beta * dE):
            lattice[x, y, z] = -lattice[x, y, z]
    return lattice

Simulation Function

Here we actually simulate the Metropolis method. First, we perform some equilibration steps to 'average out' the noise induced by the arbitrary initial state in which the system is found. After these equilibration steps, that let the system "forget" the initial condition, the configuration is statistically typical of the Boltzmann distribution for that temperature. Then, for each sweep of the algorithm, we store the magnetization and energy, which we will use to compute the following observables:

$\cdot$ Average magnetization: which we will expect, in our visual analysis, to incur in a phase transition around the critical tempoerature for the 3D Ising model.

$\cdot$ Average energy: this would have to increase proportionally to the temperature, stabilizing after the phase transition.

$\cdot$ Susceptibility $\chi=(\langle M^2\rangle-\langle M\rangle^2)/T$: this weighted magnetization variance measures how much the magnetization of the system responds to a small external magnetic field. Even when no external field is applied, it gives insights into fluctuations in magnetization and critical behavior near the phase transition.

def simulate_ising_3d(N, T, n_eq=500, n_steps=2500, J=1.0):
    rng = np.random.default_rng(SEED)
    lattice = rng.choice(np.array([-1, 1]), size=(N, N, N))
    beta = 1.0 / T

    # pre-generate all random choices for both equilibration and measurement
    total_sweeps = n_eq + n_steps
    n_trials = total_sweeps * N**3
    xs = rng.integers(0, N, size=n_trials)
    ys = rng.integers(0, N, size=n_trials)
    zs = rng.integers(0, N, size=n_trials)
    us = rng.random(n_trials)

    # Equilibration
    for sweep in range(n_eq):
        i0 = sweep * N**3
        lattice = metropolis_sweep(
            lattice, beta, J,
            xs[i0:i0+N**3], ys[i0:i0+N**3],
            zs[i0:i0+N**3], us[i0:i0+N**3]
        )

    # Observables
    mags = np.empty(n_steps)
    energies = np.empty(n_steps)
    states = []

    for sweep in range(n_steps):
        i0 = (n_eq + sweep) * N**3
        lattice = metropolis_sweep(
            lattice, beta, J,
            xs[i0:i0+N**3], ys[i0:i0+N**3],
            zs[i0:i0+N**3], us[i0:i0+N**3]
        )
        mags[sweep] = compute_magnetization(lattice)
        energies[sweep] = compute_energy(lattice, J) / lattice.size
        if sweep in [n_steps//4, n_steps//2, 3*n_steps//4]:
            states.append(lattice.copy())

    # Thermodynamic averages
    avg_M = mags.mean()
    avg_E = energies.mean()

    # Susceptibility
    chi = (np.var(mags) * N**3) / T 

    # Autocorrelation
    acf_M = autocorrelation(mags)
    
    # Distribution histograms
    hist_E, bins_E = np.histogram(energies, bins=50, density=True)
    hist_M, bins_M = np.histogram(mags, bins=50, density=True)

    return  avg_M, avg_E, chi, acf_M, (hist_E, bins_E), (hist_M, bins_M), states

Parameter Study

We now actually run the simulation for choices $N=10, 20, 30$ and a range of temperature centered at the known critical tempearature $T_C\approx 4.5$, and store the results.

Ns = [10, 20, 30]
temperatures = np.linspace(3.5, 5.0, 15)

# Run simulations in parallel over (N, T) combinations
def run_for_params(N, T):
    avg_M, avg_E, chi, acf_M, (hist_E, bins_E), (hist_M, bins_M), states = simulate_ising_3d(N, T)
    return (N, T, avg_M, avg_E, chi, acf_M, (hist_E, bins_E), (hist_M, bins_M), states)

# Flatten parameter list
tasks = [(N, T) for N in Ns for T in temperatures]

# Execute in parallel using all available cores
results_list = Parallel(n_jobs=4)(delayed(run_for_params)(N, T) for N, T in tqdm(tasks))

  0%|          | 0/45 [00:00<?, ?it/s]

# Organize results
data = {N: {'T': [], 'M': [], 'E': [], 'Chi': [], 'ACF_M': [], 'Hist_E': [], 'Hist_M': [], 'States': []} for N in Ns}
for N, T, avg_M, avg_E, chi, acf_M, (hist_E, bins_E), (hist_M, bins_M), states in results_list:
    data[N]['T'].append(T)
    data[N]['M'].append(avg_M)  
    data[N]['E'].append(avg_E)
    data[N]['Chi'].append(chi)
    data[N]['ACF_M'].append(acf_M)  
    data[N]['Hist_E'].append((hist_E, bins_E))
    data[N]['Hist_M'].append((hist_M, bins_M))
    data[N]['States'].append(states)

Plotting Results

Just to understand the temporal evolution of Metropolis-Hastings, we plot three 2D visualizations of the central cross-section slice of the lattice at different stages of the Monte Carlo run — essentially a time-lapse of how the spin configuration evolves.

fig, axes = plt.subplots(1, 3, figsize=(12, 4))
slices = [states[i][:, :, N//2] for i in [len(states)//4, len(states)//2, 3*len(states)//4]]
for ax, slc, label in zip(axes, slices, ['1/4', '1/2', '3/4']):
    ax.imshow(slc, cmap='gray')
    ax.set_title(f'State at {label} of run')
    ax.axis('off')
plt.tight_layout()
plt.show()

Then, we observe the key relation: that of magnetization as a function of temperature.

plt.figure(figsize=(8,6))

colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:purple']  
for i, N in enumerate(Ns):
    plt.scatter(data[N]['T'], data[N]['M'], color=colors[i], s=50, label=f'N={N}')

plt.axvline(4.51, color='red', linestyle='--', label='Tc ≈ 4.51')

plt.xlabel('Temperature T')
plt.ylabel('Magnetization per spin')
plt.title('3D Ising: Magnetization vs Temperature')

plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

The magnetization undergoes a sharp phase transition around the critical tempearure and this abrupt change of phase is more evident as our lattice size increases.

Now, we study the other observables.

# Susceptibility vs Temperature with peak annotation
plt.figure(figsize=(8,6))
for N in Ns:
    chi = np.array(data[N]['Chi'])
    temps = np.array(data[N]['T'])
    plt.plot(temps, chi, marker='o', label=f'N={N}')
    peak_idx = chi.argmax()
    plt.annotate(f'peak {temps[peak_idx]:.2f}', xy=(temps[peak_idx], chi[peak_idx]),
                 xytext=(temps[peak_idx]+0.05, chi[peak_idx]*1.1), arrowprops=dict(arrowstyle='->'))
plt.xlabel('Temperature T')
plt.ylabel('Susceptibility χ')
plt.title('3D Ising: Susceptibility vs Temperature')
plt.legend()
plt.grid(True)
plt.show()

At the critical temperature, the system undergoes a spontaneous symmetry-breaking transition, switching from a disordered phase (random spins) to an ordered one (aligned spins). The system is extremely sensitive to perturbations (a small field or fluctuation can tip the balance) causing huge fluctuations in magnetization, coexistence of many domains of different magnetization and spins correlated over large distances (divergent correlation length). In the last plot, this intuitive description is backed by the sharp behavior of the susceptibility around the critical temperature: when the phase transition happens, the variation is very large, otherwise, for all other temperatures, it stays very homogeneous.

plt.figure(figsize=(8, 5))

Tc = 4.51  # Critical temperature
threshold = 0.05  # How close is "near" Tc?

for T, acf in zip(data[N]['T'][::2], data[N]['ACF_M'][::2]):
    if abs(T - Tc) < threshold:
        plt.plot(acf[:100], label=f'T = {T:.2f} (≈ Tc)', color='red', linewidth=2.5, linestyle='--')
    else:
        plt.plot(acf[:100], label=f'T = {T:.2f}')

plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.title(f'Autocorrelation of Magnetization for N = {N}')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

Here we see the evolution of the autocorrelation of magnetization. When going thorugh a phase transition, from example going from solid to liquid state, a solid phase is one in which the correlations of positions are very long range, in fact they persist on arbitrary large distance: if you know the crystal structure of a solid locally, you can guess the average position of an atom even if it is very far away. Similarly here, autocorrelation decays exponentially as expected and we observe the critical slow down sensibly near to $T_c$.

Ns_to_plot = Ns  # e.g. [10, 20, 30]
Tc = 4.51
highlight_tol = 0.02  # tolerance for marking T ≈ Tc

T_indices = list(range(0, len(data[Ns_to_plot[0]]['T']), len(data[Ns_to_plot[0]]['T']) // 5))  # ~5 temps

fig, axes = plt.subplots(
    nrows=len(Ns_to_plot),
    ncols=1,
    figsize=(8, 4*len(Ns_to_plot)),
    sharex=True,
    sharey=True
)

# Collect global legend items
global_handles = []
global_labels = []
added_labels = set()

for ax, N in zip(axes, Ns_to_plot):
    Ts = data[N]['T']
    HEs = data[N]['Hist_E']
    
    for i in T_indices:
        T = Ts[i]
        hist_E, bins_E = HEs[i]
        widths = np.diff(bins_E)
        pdf = hist_E / (hist_E.sum() * widths)
        
        if abs(T - Tc) < highlight_tol:
            line = ax.step(bins_E[:-1], pdf, where='post', color='red', linewidth=2.5, label='T ≈ Tc' if 'T ≈ Tc' not in added_labels else None)
            if 'T ≈ Tc' not in added_labels:
                global_handles.append(line[0])
                global_labels.append('T ≈ Tc')
                added_labels.add('T ≈ Tc')
        else:
            label = f'T = {T:.2f}'
            line = ax.step(bins_E[:-1], pdf, where='post', label=label)
            if label not in added_labels:
                global_handles.append(line[0])
                global_labels.append(label)
                added_labels.add(label)
    
    ax.set_title(f"N = {N}")
    ax.set_ylabel("Probability density")

axes[-1].set_xlabel("Energy per spin")

# Add a single legend
fig.legend(
    global_handles,
    global_labels,
    loc='upper right',
    bbox_to_anchor=(1.18, 0.95),
    title="Temperature",
    fontsize="small"
)

fig.suptitle("Energy Histograms Across System Sizes", y=0.92)
plt.tight_layout()
plt.show()

Here what is visible is that the energy distribution is very peaked for too low/high temperatures far from the $T_c$, while, for all our choices of $N$, we notice that the energy density (red histogram) at around the transition temperature is broader. This is physically plausible, in fact around $T_c$ the systems explores a wide variety of configurations: both ordered and disordered.

plt.figure(figsize=(10, 5))
for i, N in enumerate(Ns):
    plt.plot(data[N]['T'], data[N]['M'], label=f'M, N={N}')
    plt.plot(data[N]['T'], data[N]['E'], '--', label=f'E, N={N}')
plt.xlabel('Temperature T')
plt.ylabel('Thermodynamic Averages')
plt.title('Magnetization and Energy vs Temperature for Different N')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

Finally, this plot of the thermodynamical averages of magnetization and energy confirm the global trend in their evolution and highligh their stability behavior after the phase transition: both phase portrait reach a stable plateau after passing the critical temperature.

Conclusions

Through Monte Carlo simulations of the 3D Ising model, we successfully observed the emergence of a phase transition near the critical temperature $T_c$. The magnetization, energy, and their fluctuations displayed clear signatures of this transition, including critical slowing down and broader energy distributions. Our analysis confirms that the model captures the expected physical behavior.