Physical Model and Hamiltonian: This simulation investigates the nearest-neighbor interacted spins ($S_i \in \{+1, -1\}$) located on a square lattice, under the external magnetic field. The Hamiltonian of the system is defined as:

$$ \mathcal{H} = -J \sum_{\langle i,j \rangle} S_i S_j - B_z \sum_i S_i $$

Here, $J$ is the exchange interaction constant, $\langle i,j \rangle$ represents the nearest neighbor pairs, and $B_z$ represents the external magnetic field in the z-direction. ($k_B = 1$ is assumed).

Color Coding: The spin states in the magnetic configuration visualization are color-coded as follows:

• Dark Blue: Spin Up ($S_i = +1$)
• Red: Spin Down ($S_i = -1$)

Monte Carlo (Metropolis) Update: The simulation uses the Metropolis-Hastings algorithm to reach thermodynamic equilibrium. In each trial step:

1. A random spin is selected, and the energy difference ($\Delta E$) that would result from flipping it is calculated.
2. If $\Delta E \le 0$, the move is accepted directly (the system moves to a lower or equal energy state).
3. If $\Delta E > 0$, the move is accepted with Boltzmann probability $P = e^{-\Delta E / T}$. This step allows the system to escape local minima via thermal fluctuations and preserves the principle of detailed balance.

Observables and Fluctuations:

After the system settles into thermodynamic equilibrium (following a sufficient number of Metropolis sweeps), the susceptibility ($\chi$) and specific heat ($C_v$) of the system are calculated via fluctuation theorems using the averages of magnetization $m$ and energy $e$ per spin. Let $N$ be the total number of spins ($N=100\times 100$):

• Magnetization: $$ m = \frac{1}{N} \langle \sum_i S_i \rangle $$
• Specific Heat: $$ C_v = \frac{N}{T^2} (\langle e^2 \rangle - \langle e \rangle^2) $$
• Magnetic Susceptibility: $$ \chi = \frac{N}{T} (\langle m^2 \rangle - \langle |m| \rangle^2) $$

The left chart shows the instantaneous time series of the magnetization $m$ and energy $e$ per spin (one sample every 50 MCS). The right chart shows a cumulative estimate of the fluctuation formulas above, computed from all samples collected since the last parameter change; the values settle as the number of samples grows. To keep the pre-equilibrium transient out of the averages, the first 500 MCS after each reset (thermalization) are discarded. On a finite lattice at zero field the symmetry is not spontaneously broken and $\langle m \rangle = 0$; therefore $\langle |m| \rangle$ is used in the susceptibility, as is standard.

Critical Temperature: On the infinite square lattice the phase transition occurs at $T_c = 2J/\ln(1+\sqrt{2}) \approx 2.269$ (Onsager's exact solution). The page starts at this value, and the "T = T_c" button returns to it. Near the critical point the correlation length and the autocorrelation time become very large (critical slowing down), so the χ and $C_v$ estimates can take a long time to settle. On a finite lattice the χ and $C_v$ peaks appear slightly above $T_c$ and have finite height.

Parameter Changes and Thermodynamic Memory: When the temperature (T) or magnetic field (Bz) values on the interface are modified, the simulation does not restart from a random spin configuration. Instead, the system treats its current micro-state as the initial condition and continues its evolution under the new parameters. In statistical physics, this behavior allows for the investigation of processes such as simulated annealing (gradual cooling) or thermal quench (sudden temperature drops). To prevent the mixing of data from different thermodynamic equilibrium states and the consequent corruption of variance calculations, the statistical average accumulators and charts are automatically reset upon any parameter change.