Simple Ising Model in JavaScript and D3.js

Every physicsist did at least once a simulation of the Ising model. So, why not do it twice? Let's do it in JavaScript using D3.js.

January 19, 2025

You can find ising model in almost every corner of physics. Every physicsit implemented the ising model at last once in his life. Im my case the first time I implemented the ising model was around to 2010, during my master 2. at the time I did it in C++. We are not going to do that today, we are going to try and do it in this new tool that I was planning to try for a while but I didn't had a chance too, D3.js.

The Ising model was first introduced by Wilhelm Lenz in 1920 and later solved for one dimension by Ernst Ising in 1925, and bears his name ever since. It describes the interactions of magnetic dipole moments, or “spins,” on a lattice, where each spin can take a value of +1+1 (up) or −1-1 (down). The spins are influenced by their nearest neighbors and an external magnetic field, creating a system where microscopic interactions can lead to macroscopic phenomena like phase transitions.

Up Down

The model is defined the Hamiltonian, which represents the system's total energy and spell how each individual spin interect with its neighbours. For a lattice with spins \sigma_i, the Hamiltonian thake the expression of

H=−J∑⟨i,j⟩σiσj−h∑iσi.H = -J \sum_{\langle i, j \rangle} \sigma_i \sigma_j - h \sum_i \sigma_i.

Here, JJ is the coupling constant dictating the interaction strength between spins: positive JJ encourages alignment (ferromagnetic behavior), while negative JJ favors anti-alignment (antiferromagnetic behavior). The term h represents the external magnetic field, and ⟨i,j⟩\langle i, j \rangle signifies summation over neighboring pairs. The competition between these terms and the thermal fluctuations at a given temperature determines the system's behavior.

A bit more Ising theory

This section fills in the exact 2D results and the simulation details you can validate against. We consider an L×LL\times L square lattice with spins σi∈{±1}\sigma_i\in\{\pm1\}, nearest-neighbour coupling JJ, and (optional) external field hh. The Hamiltonian is

H({σ})=−J ⁣∑⟨i,j⟩σiσj−h ⁣∑iσi.H(\{\sigma\}) = -J\!\sum_{\langle i,j\rangle} \sigma_i\sigma_j - h\!\sum_i \sigma_i.

It is often convenient to work with dimensionless couplings

K≡βJ,h~≡βh,β=(kBT)−1.K \equiv \beta J,\qquad \tilde h \equiv \beta h,\qquad \beta = (k_B T)^{-1}.

On a square lattice the coordination number is z=4z=4. For simulations we usually impose periodic boundary conditions (PBCs). What happens in 1D vs 2D?

  • 1D Ising (h=0h=0): no finite-temperature phase transition; Tc=0T_c=0.
  • 2D square-lattice Ising (h=0h=0) is exactly solvable and exhibits a finite TcT_c.

Exact results in 2D (h=0h=0)

Critical temperature. Onsager's solution gives the critical coupling via

sinh⁡(2Kc)=1  ⇒  Kc=12ln⁡(1+2)≈0.44068679…\sinh(2K_c)=1 \;\Rightarrow\; K_c=\tfrac12\ln(1+\sqrt2)\approx 0.44068679\ldots

If you set J=kB=1J=k_B=1, this corresponds to

Tc=1Kc≈2.2691853 .T_c = \frac{1}{K_c} \approx 2.2691853\,.

Spontaneous magnetization for T<TcT<T_c. Yang showed that the exact magnetization is

m(T)=[1−sinh⁡−4(2K)]1/8,K=J/(kBT),T<Tc.m(T) = \left[1-\sinh^{-4}(2K)\right]^{1/8},\quad K=J/(k_B T),\quad T<T_c.

Near TcT_c, m∼(Tc−T)βm\sim (T_c-T)^{\beta} with critical exponent β=1/8\beta=1/8.

Other famous 2D critical exponents.

α=0  (log divergence of cV),γ=74,ν=1,η=14.\alpha=0\;\text{(log divergence of $c_V$)},\quad \gamma=\tfrac{7}{4},\quad \nu=1,\quad \eta=\tfrac{1}{4}.

These numbers are handy for checking finite‑size scaling.

Observables to measure

With N=L2N=L^2:

m=1N∑iσi,χ=βN (⟨m2⟩−⟨m⟩2),m = \frac{1}{N}\sum_i\sigma_i,\qquad \chi = \beta N\,(\langle m^2\rangle-\langle m\rangle^2), ϵ=HN,cV=β2N(⟨H2⟩−⟨H⟩2).\epsilon = \frac{H}{N},\qquad c_V = \frac{\beta^2}{N}\big(\langle H^2\rangle-\langle H\rangle^2\big).

The Binder cumulant is a very useful dimensionless quantity for locating TcT_c with finite lattices:

U4=1−⟨m4⟩3⟨m2⟩2.U_4 = 1-\frac{\langle m^4\rangle}{3\langle m^2\rangle^2}.

For different LL, the U4(T)U_4(T) curves cross near TcT_c.

Metropolis update — explicit ΔE\Delta E

Flipping a single spin σi→−σi\sigma_i\to-\sigma_i changes the energy by

ΔE=2J σi ⁣∑j∈nn(i) ⁣σj+2h σi.\Delta E = 2 J\,\sigma_i\!\sum_{j\in \text{nn}(i)}\!\sigma_j + 2 h\,\sigma_i.

Accept the flip with probability min⁡{1, e−βΔE}\min\{1,\,e^{-\beta\Delta E}\}. On a square lattice with h=0h=0, the local field can only take values in {−4J,−2J,0,2J,4J}\{-4J,-2J,0,2J,4J\}, so ΔE∈{−8J,−4J,0,4J,8J}\Delta E\in\{-8J,-4J,0,4J,8J\}. You can precompute the five Boltzmann factors to avoid repeated exp calls.

One Monte Carlo sweep (MCS): attempt N=L2N=L^2 single‑spin flips. Discard an initial thermalization window (e.g., 102 ⁣− ⁣10310^2\!-\!10^3 MCS), then sample every few sweeps to reduce autocorrelation.

Finite‑size scaling checklist

  • Pseudocritical temperature shifts as
Tc(L)≈Tc+a L−1/ν(ν=1  in 2D).T_c(L)\approx T_c + a\,L^{-1/\nu}\quad (\nu=1\;\text{in 2D}).
  • Peaks of χ\chi grow χmax⁡∝Lγ/ν=L7/4\chi_{\max}\propto L^{\gamma/\nu}=L^{7/4}.
  • Specific‑heat peak grows roughly like log⁡L\log L (since α=0\alpha=0).

Practical validation targets (with J=kB=1J=k_B=1)

  • Critical coupling: Kc≈0.44068679K_c\approx 0.44068679. Equivalent statement often used in code: βeff,c=Jβc=Kc\beta_{\text{eff},c}=J\beta_c=K_c.
  • Critical temperature: Tc≈2.2691853T_c\approx 2.2691853.
  • Magnetization check: pick T=2.0T=2.0 (i.e., K=0.5K=0.5) and compare your measured mm to the exact [1−sinh⁡−4(1)]1/8≈0.953[1-\sinh^{-4}(1)]^{1/8}\approx 0.953 as a sanity check.

Understanding the Ising model requires calculating quantities such as the total energy, magnetization, and correlation functions. In principle, these can be derived from the partition function

Z=∑{σ}e−βHZ = \sum_{\{\sigma\}} e^{-\beta H}

where β=1/kBT\beta = 1/k_B T is the inverse temperature. The partition function sums over all possible spin configurations, with each weighted by its Boltzmann factor e−βHe^{-\beta H}. From ZZ, macroscopic quantities like free energy, specific heat, and magnetization can be computed. However, this summation scales exponentially with the number of spins, making exact solutions infeasible for large systems.

To address this computational challenge, Monte Carlo methods are employed. These numerical techniques allow us to sample the most probable spin configurations rather than exhaustively enumerating all possibilities. A widely used algorithm for simulating the Ising model is the Metropolis algorithm, which mimics the thermal fluctuations of a physical system. The algorithm begins by initializing a lattice of spins, typically in a random configuration. At each step:

  1. a spin is chosen at random and flipped, changing its value from +1+1 to −1-1 or vice versa.
  2. The energy difference \Delta E resulting from this flip is computed. If the flip reduces the energy (ΔE≤0)(\Delta E \leq 0), it is accepted. Otherwise, it is accepted with a probability proportional to e−βΔEe^{-\beta \Delta E}, allowing the system to escape local energy minima and explore a range of configurations.

The classical Ising model assumes spins are discrete variables, taking values of +1+1 or −1-1. This contrasts with the quantum Ising model, where spins are treated as quantum mechanical operators capable of existing in superpositions. The quantum model introduces additional complexity, such as a transverse field term, which drives quantum fluctuations and leads to phenomena like quantum phase transitions. While the classical model can be simulated with Monte Carlo techniques, the quantum version often requires advanced methods like Quantum Monte Carlo or tensor network approaches.

Simulating the Ising model provides insights into critical phenomena and phase transitions. By adjusting parameters like temperature, coupling strength, and external field, one can observe how spins transition from an ordered state, where most are aligned, to a disordered state dominated by thermal fluctuations. The transition point, or critical temperature, marks the onset of long-range correlations, a hallmark of second-order phase transitions.

Implementing the Ising model computationally involves defining the lattice, initializing spins, and simulating their evolution using Monte Carlo updates. By periodically measuring quantities like magnetization and energy, one can study the system's response to varying parameters. Visualization tools, such as D3.js, can enhance this analysis by dynamically displaying the lattice and its evolution over time.

The Ising model's simplicity belies its depth, making it a valuable tool for exploring fundamental principles of statistical mechanics. Its applications extend far beyond magnetism, offering insights into problems as diverse as neural network dynamics, social behavior, and even protein folding. Through simulations, we not only understand the theoretical framework but also uncover the model's rich phenomenology, bridging the gap between microscopic rules and macroscopic order.

Building the sumulation class in JavaScript

We start by defining class that will host our simulation, let call it IsingSimulation. This class will contain the parameters of the model, such as the grid size, coupling constant, external field, and temperature. It will also store the spin configuration and provide methods for calculating energy and running Monte Carlo updates.

class IsingSimulation {
  constructor(N, J, h, temperature) {
    this.N = N; // Grid size (NxN)
    this.Nsites = N * N; // Total number of sites
    this.J = J; // Coupling constant
    this.h = h; // External magnetic field
    this.temperature = temperature; // Temperature of the system
    this.spins = new Array(this.Nsites)
      .fill(0)
      .map(() => (Math.random() > 0.5 ? 1 : -1)); // Initialize spins randomly
  }
}

You will notice that the spins are initialized randomly, with a 50% chance of being up or down. This randomness ensures the system explores a wide range of configurations during the simulation.

Now it will important to define a function within this class that will calculate the energy of the system. The energy of the Ising model is given by the Hamiltonian, which sums over all pairs of neighboring spins:

calculateEnergy() {
  let energy = 0;
  for (let i = 0; i < this.Nsites; i++) {
    const x = i % this.N;
    const y = Math.floor(i / this.N);
    let neighborSum = 0;
 
    // Loop over neighbors
    for (let dx = -1; dx <= 1; dx++) {
      for (let dy = -1; dy <= 1; dy++) {
        if (dx === 0 && dy === 0) continue; // Skip self
        const nx = (x + dx + this.N) % this.N;
        const ny = (y + dy + this.N) % this.N;
        const j = nx + ny * this.N;
        neighborSum += this.spins[j];
      }
    }
    energy -= this.J * this.spins[i] * neighborSum - this.h * this.spins[i];
  }
  return energy / 2; // To avoid double-counting
}

In this function, neighbors are identified using modular arithmetic to implement periodic boundary conditions. For each site, interactions with neighbors and the external field contribute to the energy.

Now the hart of the simulation, the Monte Carlo method, using the Metropolis algorithm, evolves the system:

runMonteCarlo(numSteps) {
  let energy = this.calculateEnergy();
  for (let step = 0; step < numSteps; step++) {
    const i = Math.floor(Math.random() * this.Nsites); // Choose a random site
    const newSpin = -this.spins[i]; // Flip the spin
    this.spins[i] = newSpin;
 
    const newEnergy = this.calculateEnergy();
    const deltaE = newEnergy - energy;
 
    // Accept or reject the flip based on energy change
    if (deltaE > 0 && Math.random() > Math.exp(-deltaE / this.temperature)) {
      this.spins[i] = -newSpin; // Reject the flip
    } else {
      energy = newEnergy; // Accept the flip
    }
  }
}

Key Points to nitice here. A random spin is flipped. The energy change, ΔE\Delta E, determines if the flip is accepted. Flips that reduce energy are always accepted; others are probabilistically accepted based on the Boltzmann factor. This stochastic process mimics thermal fluctuations in the physical system. The number of Monte Carlo steps determines the simulation's duration and the extent of exploration in phase space.

Visualizing the Spins with D3.js

Using D3.js, the spins are displayed as a grid of colored circles:

const N = 20; // Grid size
const width = 200; // Visualization dimensions
const height = 200;
const cellSize = width / N;
 
const svg = d3
  .select("#isingModel")
  .append("svg")
  .attr("width", width)
  .attr("height", width);
 
function visualizeSpins() {
  const spins = simulation.spins;
  const cells = svg
    .selectAll("circle")
    .data(spins)
    .join("circle")
    .attr("cx", (d, i) => (0.5 + (i % N)) * cellSize)
    .attr("cy", (d, i) => (0.5 + Math.floor(i / N)) * cellSize)
    .attr("r", cellSize * 0.4)
    .attr("fill", (d) => (d === 1 ? "steelblue" : "orange"));
 
  cells.exit().remove();
}

Spins are represented by circles: +1+1 red and −1-1 is a transparent white. The visualization updates dynamically as the simulation progresses. Positions are calculated based on the grid index.

Finally, once can change the dymanically the parameters of the ising model using three sliders allow users to control temperature TT, coupling strength JJ, and external field hh:

function updateTemperature(value) {
  simulation.temperature = parseFloat(value);
  document.getElementById("tValue").textContent = value;
}

Similar functions adjust JJ and hh. The interface dynamically updates the simulation in real-time.

Conclusion

The Ising model is a fascinating way to explore phase transitions and critical phenomena. By implementing the model in JavaScript and embedding in my Next.js website, the idea is to make these concepts more tangible and interactive. One could whish to extend the simulation by adding features like heat maps, energy plots, or magnetization curves. The goal is to provide an intuitive way to engage with complex theoretical concepts and appreciate the beauty of statistical mechanics. I will try to do this in the future if time allows.

References

[^ref]