Physical Model and Hamiltonian: Including 3D external magnetic fields and chiral interactions, the full Hamiltonian for continuous spin vectors ($|\mathbf{S}_i| = 1$) is defined as:
$$ \mathcal{H} = -J \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j - \sum_i \mathbf{B} \cdot \mathbf{S}_i - \sum_{\langle i,j \rangle} \mathbf{D}_{ij} \cdot (\mathbf{S}_i \times \mathbf{S}_j) $$
Dzyaloshinskii-Moriya Interaction (DMI): The antisymmetric exchange stabilizes topological textures. The orientation of the DMI vector $\mathbf{D}_{ij}$ depends on the structural inversion asymmetry:
Topological Charge ($Q$): The topological invariant $Q$ counts how many times the spin configuration wraps the unit sphere. In the continuum it is $Q = \frac{1}{4\pi}\int \mathbf{S}\cdot(\partial_x \mathbf{S} \times \partial_y \mathbf{S})\,dx\,dy$. On the discrete lattice this simulation uses the Berg–Lüscher definition: each square plaquette is split into two triangles, and the solid angle spanned by the three spins of each triangle on the unit sphere is computed:
$$ \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 boundary conditions the lattice is a torus, and this definition gives an integer for every configuration, independent of temperature or of how sharp the gradients are. The $Q$ chart therefore changes in steps. Each step is an event in which a skyrmion is created or destroyed on the lattice. This process is topologically forbidden in the continuum but possible at the lattice scale, in a single cell where neighbouring spins are nearly antiparallel. The simple triple-product approximation $\frac{1}{4\pi}\sum \mathbf{S}_i\cdot(\mathbf{S}_{i+\hat x}\times\mathbf{S}_{i+\hat y})$ is accurate only for slowly varying textures; for skyrmions a few lattice constants across it clearly undercounts.
Sign and Length Scale: For $B_z > 0$ the background points along $+z$, and each skyrmion with a $-z$ core contributes $Q = -1$; thus $|Q|$ is the number of skyrmions. The characteristic helix period is $\lambda \approx 2\pi J/|D|$ for small $D$, and on the square lattice it is exactly given by $\tan(2\pi a/\lambda) = |D|/J$ ($\lambda = 8a$ for $D = J$). When $\lambda$ does not fit the system size ($L=30$), the periodic boundaries compress or stretch the texture. This distortion affects the shape and number of skyrmions, but it does not make $Q$ fractional.
Monte Carlo Algorithm: The simulation uses the Metropolis–Hastings algorithm. A trial spin is generated inside a "cone" around the current spin as $\mathbf{S}' = (\mathbf{S} + \delta\,\mathbf{r})/|\mathbf{S} + \delta\,\mathbf{r}|$. Here $\mathbf{r}$ is uniformly distributed on the unit sphere and is obtained with Marsaglia's method (1972): $u, v \in [-1, 1]$ are drawn until $s = u^2 + v^2 < 1$, and $\mathbf{r} = (2u\sqrt{1-s},\, 2v\sqrt{1-s},\, 1-2s)$. This proposal is symmetric, so detailed balance is preserved. The cone width $\delta$ is slowly tuned after every sweep to keep the acceptance rate near 50%. Fully random trial spins are almost never accepted at low temperature (about 3% at $T = 0.15$); cone moves prevent the system from artificially freezing.
Update Frequency (MCS/frame): Time is measured in Monte Carlo Sweeps (MCS). This parameter sets how many MCS are performed between two screen redraws. A higher value does not reach equilibrium in fewer MCS, but it advances the system faster in real (wall-clock) time and makes successive visible frames less correlated.
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, enabling the direct simulation of thermal quenches or simulated annealing processes. However, to prevent statistical contamination between different equilibrium states, the magnetization arrays, charge arrays, and visualization charts are automatically cleared upon any parameter alteration.
Suggested Experiments: The page starts from a random state with $T = 0.15$, $B_z = 0.5$, $D = 1$ (Néel). Within a few hundred MCS a texture of skyrmions with $-z$ cores forms, with a topological charge of about $|Q| \approx 15$–$20$; the "Skyrmion example" button returns to this starting point. Setting $B_z = 0$ and pressing "Random start" produces labyrinth-like spiral domains instead of skyrmions, and $Q$ stays close to zero. Existing skyrmions, by contrast, survive for a long time even after the field is switched off, a direct illustration of topological protection. As $B_z$ increases the skyrmions shrink and their number decreases; for $B_z \gtrsim 1$ the system becomes ferromagnetic. Note that for the Néel type the in-plane arrows point radially with respect to the skyrmion centre, while for the Bloch type they are tangential.
Visualization and Color Coding: 3D spin vectors are projected onto the 2D lattice plane. The background color of a cell represents the $S_z$ component (Positive $\to$ Dark Blue, Negative $\to$ Red). The black arrow inside represents the $(S_x, S_y)$ projection.