Physical Model and Hamiltonian: With a 3D external magnetic field $\mathbf{B} = (B_x, B_y, B_z)$ and continuous spin vectors $|\mathbf{S}_i| = 1$, the Hamiltonian of the system is:

$$ \mathcal{H} = -J \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j - \sum_i \mathbf{B} \cdot \mathbf{S}_i $$

Monte Carlo Update Algorithm and Marsaglia Sampling: The simulation uses the Metropolis–Hastings algorithm. Trial spins are generated inside a "cone" around the current spin. First, a random vector $\mathbf{r}$ uniformly distributed on the unit sphere is generated with Marsaglia's method (1972), which avoids trigonometric functions: $u, v \in [-1, 1]$ are drawn until $s = u^2 + v^2 < 1$. The trial spin is then:

$$ \mathbf{r} = \left(2u\sqrt{1-s},\; 2v\sqrt{1-s},\; 1-2s\right), \qquad \mathbf{S}' = \frac{\mathbf{S} + \delta\,\mathbf{r}}{|\mathbf{S} + \delta\,\mathbf{r}|} $$

This proposal is symmetric (its probability depends only on $\mathbf{S}\cdot\mathbf{S}'$), so detailed balance is preserved. The cone width $\delta$ is slowly tuned after every sweep to keep the acceptance rate near 50%; the current values are shown in the control panel. (Because the tuning happens only between sweeps, it introduces a negligible bias for this demonstration.) Fully random trial spins on the sphere are almost never accepted at low temperature and the system artificially "freezes"; cone moves prevent this.

Visualization and Color Coding: The lattice projection maps 3D spin vectors onto a 2D plane:

• Z-Component ($S_z$): Represented by the background color of the cell. Positive values ($S_z \to +1$) transition to Dark Blue, while negative values ($S_z \to -1$) transition to Red. Intermediate values near the $x-y$ plane approach white.
• X-Y Plane Projection: Represented by the black arrow inside each cell. The orientation and scalar length of the arrow directly depict the $(S_x, S_y)$ vector components.

Thermodynamic Memory: Modifying the temperature or magnetic field parameters during the simulation does not reset the lattice to a random configuration. The system retains its current micro-state and resumes thermodynamic evolution under the new conditions. This feature enables the direct simulation of thermal quenches or simulated annealing processes. However, to prevent statistical contamination between different equilibrium states, the magnetization arrays and visualization charts are automatically cleared upon any parameter alteration.

Update Frequency (MCS/frame): In Monte Carlo simulations, time is measured in Monte Carlo Sweeps (MCS), where one MCS corresponds to $N$ trial spin updates. This parameter sets how many MCS are performed between two screen redraws. A higher value does not make the system reach equilibrium in fewer MCS. Because the rendering cost is spread over fewer frames, however, the system advances faster in real (wall-clock) time and successive visible frames are less correlated. A value of 1 shows closely neighbouring micro-states of the Markov chain one by one.

Mermin-Wagner Theorem and Finite-Size Effects: For an isotropic 2D Heisenberg model with short-range interactions at zero external field ($\mathbf{B} = 0$), the Mermin-Wagner theorem rigorously dictates that continuous rotational symmetry cannot be spontaneously broken at any finite temperature ($T > 0$). Consequently, in the thermodynamic limit ($N \to \infty$), the system cannot sustain long-range ferromagnetic order, and the expectation value of global magnetization is zero. Any apparent macroscopic ordering or non-zero magnetization observed in this specific simulation at low temperatures is strictly a finite-size effect. It occurs because the magnetic correlation length exceeds the limited spatial dimensions of the simulated lattice ($L=30$), trapping the system in a pseudo-ordered state.

Observables: The line chart tracks the time evolution (in Monte Carlo Sweeps) of the global magnetization components:

$$ M_{\alpha} = \frac{1}{N} \sum_i S_i^{\alpha}, \quad \alpha \in \{x, y, z\} $$