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 employs the Metropolis-Hastings algorithm. In each trial step, a random continuous spin must be generated uniformly on the unit sphere. To achieve high computational efficiency without relying on expensive trigonometric functions, the algorithm utilizes Marsaglia's rejection method (1972). Two random variables $u, v \in [-1, 1]$ are drawn until $s = u^2 + v^2 < 1$. The trial spin components are then evaluated as:

$$ S_x = 2u\sqrt{1-s}, \quad S_y = 2v\sqrt{1-s}, \quad S_z = 1-2s $$

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 strictly measured in Monte Carlo Sweeps (MCS), where one MCS corresponds to $N$ independent trial spin updates. The browser's visualization process (rendering thousands of vectors) creates a severe computational bottleneck. This parameter explicitly decouples the physical Markov chain evolution from the visual rendering rate. A higher value (e.g., 50) rapidly pushes the system through phase space, effectively minimizing the autocorrelation between successive visual frames and accelerating convergence to thermodynamic equilibrium. Conversely, a value of 1 visualizes highly correlated, adjacent micro-states in the Markov chain.

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\} $$