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:

• Néel Type (Interfacial): Stabilizes chiral Néel skyrmions or spin cycloids. Defined as $\mathbf{D}_{ij} = D (\hat{u}_{ij} \times \hat{z})$.
• Bloch Type (Bulk): Stabilizes Bloch skyrmions or spin helices. Defined as $\mathbf{D}_{ij} = D \hat{u}_{ij}$.

Topological Charge ($Q$) and Fractional Values: The topological invariant $Q$ counts how many times the spin configuration wraps around the unit sphere. It is numerically approximated on the discrete lattice using the triple scalar product of neighboring spins:

$$ Q \approx \frac{1}{4\pi} \sum_{i,j} \mathbf{S}_{i,j} \cdot (\mathbf{S}_{i+1,j} \times \mathbf{S}_{i,j+1}) $$

In a continuous thermodynamic limit at absolute zero, $Q$ must be a strictly quantized integer ($0, \pm 1, \pm 2, \dots$). However, in this simulation, the topological charge time-series may frequently display fractional (non-integer) values. This artifact stems from two fundamental reasons:

1. Commensurability and PBC Truncation: If the intrinsic skyrmion pitch length ($\lambda \approx 2\pi J / |D|$) is not commensurate with the finite lattice size ($L=30$), the Periodic Boundary Conditions (PBC) artificially slice the skyrmion at the boundaries. An incomplete, frustrated topological texture mathematically evaluates to a fractional $Q$.
2. Discrete Finite-Difference Approximation: The triple scalar product is only an approximation of the continuous spatial integral. When the DMI is strong or thermal fluctuations ($T > 0$) are high, the angle between adjacent spins becomes large. In such regimes of sharp spin gradients, this approximation breaks down. A rigorous, computationally expensive geometric method (e.g., Berg-Lüscher) would be required to recover exact integers.

Monte Carlo Algorithm and Marsaglia Sampling: The simulation utilizes the Metropolis-Hastings algorithm. To generate uniform random trial spins on the unit sphere without the heavy computational cost of trigonometric functions, Marsaglia's rejection method (1972) is employed. Two variables $u, v \in [-1, 1]$ are drawn until $s = u^2 + v^2 < 1$. The new spin components are exactly mapped as $S_x = 2u\sqrt{1-s}, S_y = 2v\sqrt{1-s}, S_z = 1-2s$.

Update Frequency (MCS/frame): Time is strictly measured in Monte Carlo Sweeps (MCS). To bypass the browser's rendering bottleneck, this parameter decouples the physical Markov chain from visual frame rates. A higher value accelerates the system through phase space by performing multiple MCS invisibly, minimizing autocorrelation between visible frames and speeding up thermodynamic convergence.

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.

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.