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{tot}} - \alpha \, \mathbf{S}_i \times (\mathbf{S}_i \times \mathbf{h}_i^{\text{tot}}) $$
Here $\mathbf{h}_i^{\text{tot}} = \mathbf{h}_i^{\text{eff}} + \mathbf{h}_i^{\text{th}}$ is the sum of the effective and thermal fields.
Thermal Noise and Langevin Dynamics: Finite temperature ($T > 0$) is modelled by adding a stochastic thermal field ($\mathbf{h}_i^{\text{th}}$) to the effective field. This white noise is uncorrelated in space and time, and its strength is set by the fluctuation–dissipation theorem. The Gilbert-form correlation
$$\langle h^{\text{th}}_\mu(t)\, h^{\text{th}}_\nu(t')\rangle = (2\alpha k_B T/\gamma\mu_s)\,\delta_{\mu\nu}\,\delta(t-t')$$
acquires a factor $1/(1+\alpha^2)$ in units of $\tau$. The standard deviation of the Gaussian noise added at each step $d\tau$ is therefore:
$$ \sigma = \sqrt{\frac{2\alpha T}{(1+\alpha^2)\, d\tau}} $$
Gaussian random numbers are generated from uniform random numbers using the Box–Muller transform. (Check: for non-interacting spins this calibration makes $\langle S_z \rangle$ follow the Langevin function $\coth(B/T) - T/B$.)
Stratonovich Integration: The thermal field enters the equation multiplied by the spin vector (multiplicative noise). The interpretation of the stochastic integral therefore matters: physically, the LLG equation must be read in the Stratonovich sense. The simple Euler scheme corresponds to the Itô interpretation and, without an additional drift correction, leads to the wrong equilibrium distribution. This simulation uses the Heun (predictor–corrector) method with the thermal field held fixed during both stages; the Heun method converges to the Stratonovich solution. The spins are renormalized after each stage to preserve $|\mathbf{S}_i| = 1$. The finite step $d\tau = 0.02$ introduces a small systematic error; the Boltzmann distribution is recovered exactly in the limit $d\tau \to 0$.
$C_6$ Geometry and DMI Vectors: The simulation is built on a triangular lattice in which each spin has 6 neighbours arranged at $\pi/3$ angles. The DMI vectors ($\mathbf{D}_{ij}$) are defined through the unit vectors pointing to the neighbours ($\hat{u}_{ij}$): $\mathbf{D}_{ij} = D(\hat{u}_{ij} \times \hat{z})$ for Néel type and $\mathbf{D}_{ij} = D\hat{u}_{ij}$ for Bloch type. Since the Hamiltonian is
$$\mathcal{H} = -J\sum_{\langle i,j\rangle}\mathbf{S}_i\cdot\mathbf{S}_j - \sum_{\langle i,j\rangle}\mathbf{D}_{ij}\cdot(\mathbf{S}_i\times\mathbf{S}_j) - \sum_i \mathbf{B}\cdot\mathbf{S}_i$$
the effective field acting on each spin is:
$$ \mathbf{h}_i^{\text{eff}} = -\frac{\partial \mathcal{H}}{\partial \mathbf{S}_i} = \sum_{j} \left( J\,\mathbf{S}_j + \mathbf{S}_j \times \mathbf{D}_{ij} \right) + \mathbf{B} $$
System Size and Periodic Boundaries: The simulation runs on an $L \times L$ ($L=30$) rhombus-shaped lattice with periodic boundary conditions (PBC). The characteristic helix period (in lattice constants) is of order $\lambda \sim 2\pi J/|D|$. When $\lambda$ does not fit the system size, the periodic boundaries compress or stretch the texture, which affects the shape, number and arrangement of the skyrmions. Well-defined skyrmions require $\lambda$ to be sufficiently smaller than $L$.
Topological Charge ($Q$): The topological charge is computed with the Berg–Lüscher definition over the triangles that completely tile the $C_6$ lattice (one up and one down triangle per lattice site). The solid angle spanned by the three spins of each triangle on the unit sphere is:
$$ \Omega_{\triangle} = 2\arctan\frac{\mathbf{S}_1\cdot(\mathbf{S}_2\times\mathbf{S}_3)}{1+\mathbf{S}_1\cdot\mathbf{S}_2+\mathbf{S}_2\cdot\mathbf{S}_3+\mathbf{S}_3\cdot\mathbf{S}_1}, \qquad Q = \frac{1}{4\pi}\sum_{\triangle}\Omega_{\triangle} $$
With periodic boundaries this definition gives an integer for every configuration; the $Q$ chart steps only when a skyrmion is created or destroyed at the lattice scale. For $B_z > 0$ each skyrmion with a $-z$ core contributes $Q = -1$. (Note: at small angles the plain triple product $\mathbf{S}_1\cdot(\mathbf{S}_2\times\mathbf{S}_3)$ is about twice the solid angle, so summing it over all triangles overestimates $Q$ by roughly a factor of two.)
Time Step and Thermodynamic Memory: The "Integrations/Frame" parameter sets how many Heun steps of size $d\tau = 0.02$ are performed between two screen redraws; the horizontal axis of the charts shows the dimensionless time $\tau$. When a parameter is changed, the instantaneous spin state is kept and the evolution continues from where it left off (so quench or annealing experiments are possible); only the time-series charts are reset.
Suggested Experiment: The page starts from a random state with $T = 0.1$, $B_z = 1$, $D = 1$ (Néel), and within a few tens of $\tau$ skyrmions with $-z$ cores form in a ferromagnetic background (the "Skyrmion example" button returns to this starting point). Starting instead from the uniform state with "Ferromagnetic (+z) start", no skyrmions form at the same parameters: the uniform state is metastable, and small thermal fluctuations cannot overcome the topological barrier. This comparison directly illustrates the topological protection of skyrmions.
Visualization and Color Coding: The lattice projection maps 3D spin vectors onto the 2D plane: