Dimensionlessization of the LLG Equation: To enhance the numerical stability of the simulation, physical time is dimensionlessized as $\tau = t / t_0$ using the unit $t_0 = (1+\alpha^2)/(\gamma J)$. Through energy scaling ($J=1$), the effective field derived from the Hamiltonian ($\mathbf{h}_i^{\text{eff}} = -\partial \mathcal{H} / \partial \mathbf{S}_i$) is transformed into a dimensionless form free of micromagnetic constants. The equation containing deterministic precession and damping terms takes the following form:
$$ \frac{d\mathbf{S}_i}{d\tau} = -\mathbf{S}_i \times \mathbf{h}_i^{\text{toplam}} - \alpha \, \mathbf{S}_i \times (\mathbf{S}_i \times \mathbf{h}_i^{\text{toplam}}) $$
Thermal Noise and Langevin Dynamics: Finite temperature ($T > 0$) is modeled by adding a stochastic thermal noise field ($\mathbf{h}_i^{\text{th}}$) to the effective field, in accordance with the Fluctuation-Dissipation theorem. This white noise has zero correlation in space and time. In the integration performed via time discretization, the standard deviation of the random Gaussian noise added at each step is determined by $\sigma = \sqrt{2 \alpha T / d\tau}$. To overcome the limitations of browser-based pseudo-random number generators, normal (Gaussian) distributed 3D vectors are generated using the Box-Muller transform.
Stratonovich Integration: The thermal field appears multiplicatively with the spin vector in the equation. When the Euler method (Itô limit) is used, the unit spin constraint of the system ($|\mathbf{S}|=1$) is violated. Stratonovich calculus is required to converge to the physical Boltzmann distribution. The algorithm utilizes the Heun (Predictor-Corrector) integrator and guarantees thermodynamic equilibrium by keeping the thermal noise vector ($\mathbf{h}_i^{\text{th}}$) constant during both integration stages.
$C_6$ Geometry and DMI Vectors: The simulation is built on a triangular lattice where each spin has 6 neighbors arranged at $\pi/3$ angles. DMI vectors ($\mathbf{D}_{ij}$) are analytically projected via the unit vectors of the neighbors ($\hat{u}_{ij}$). It follows the orientations $\mathbf{D}_{ij} = D(\hat{u}_{ij} \times \hat{z})$ for Néel-type DMI and $\mathbf{D}_{ij} = D\hat{u}_{ij}$ for Bloch-type DMI.
System Size, PBC, and Topological Frustration: The simulation operates on an $L \times L$ sized ($L=30$) lattice tiled in a rhombus form. Edge effects are eliminated by using Periodic Boundary Conditions (PBC). However, when the intrinsic skyrmion pitch length ($\lambda \approx 2\pi J / |D|$) does not fit the finite lattice size, PBC artificially truncates the skyrmion, causing topological frustration. To observe distinct skyrmionic phases, the $|D|$ value must satisfy the condition $\lambda \le L$.
Topological Charge ($Q$) and Finite Difference Deviation: The topological invariant is calculated in discrete form by summing the scalar triple products of the triangular sub-cells completely tiling the $C_6$ lattice:
$$ Q \approx \frac{1}{4\pi} \sum_{\triangle} \mathbf{S}_1 \cdot (\mathbf{S}_2 \times \mathbf{S}_3) $$
In continuous space, $Q$ is an absolutely quantized integer. However, the scalar triple product in this simulation is treated via a finite difference approximation. In high-temperature fluctuations or sharp gradient regimes where DMI dramatically bends the angle between field vectors, this approach becomes inadequate, and fractional deviations are observed in the $Q$ chart.
Time Steps (dt) and Thermodynamic Memory: To alleviate the client-side processing load, physical time ($\tau$) and visualization frames are decoupled. The "Integrations/Frame" parameter determines how many Heun integrations of size $d\tau = 0.02$ are performed in the background before drawing to the screen. During parameter changes, the instantaneous spin micro-state of the system is preserved, and thermodynamic evolution continues from where it left off (simulated annealing); however, to maintain simplicity in the graphical representation, the statistical pools are reset, meaning the time series charts restart when parameters are changed.
Visualization and Color Coding: The lattice projection maps 3D spin vectors onto the 2D plane: