Structure-Preserving Spin-Lattice Integrator

• Structure-Preserving Spin-Lattice Integrator is a numerical scheme that maintains key invariants such as spin norm, momentum, and energy in coupled spin-lattice dynamics.
• It employs Hamiltonian splitting, exact rotational updates, and conservative thermostats to ensure long-term accuracy and stability in simulations.
• The method underpins large-scale simulations and machine-learning potentials, enabling precise study of magnetic phase transitions and spin-wave behavior.

A structure-preserving spin-lattice integrator is a numerical scheme for spin systems or coupled spin-lattice dynamics whose discrete update is constructed to retain the geometric properties of the underlying equations of motion. In recent work, this designation covers Hamiltonian splittings for positions, momenta, and classical spins; exact rotation maps for Lie-Poisson spin precession; conservative stochastic thermostats that preserve a conserved order parameter; and norm-preserving implicit schemes for Landau-Lifshitz evolution. Depending on the formulation, the preserved structures include the spin norm $|\mathbf S_i|=1$, global magnetization, total linear and angular momentum, time-reversibility, symplecticity, bounded long-time energy error, or exact sampling of the Gibbs ensemble $P(\sigma)\propto e^{-H/T}$ (Tranchida et al., 2018, Weißenhofer et al., 2022, Cavagna et al., 2022, Huang et al., 15 Jun 2025, Chen et al., 12 Jun 2026, Xie et al., 11 Feb 2026).

1. Hamiltonian formulations and geometric setting

The standard spin-lattice formulation couples atomic positions $\mathbf r_i$, momenta $\mathbf p_i$, and classical spins $\mathbf s_i$ or $\mathbf S_i$. A representative Hamiltonian is

$$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$

with equations of motion

$$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$

In this setting, the lattice sector is canonical, whereas the spin sector is non-canonical and of Lie-Poisson type (Tranchida et al., 2018, Chen et al., 12 Jun 2026).

A more explicit coupled Hamiltonian combines a magnetic term $\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$, a kinetic term, and a mechanical potential $V(r_{ij})$, with the magneto-mechanical coupling entering through the distance dependence of the exchange $P(\sigma)\propto e^{-H/T}$0. The resulting equations contain the mixed canonical/non-canonical structure

$P(\sigma)\propto e^{-H/T}$1

with $P(\sigma)\propto e^{-H/T}$2 (Tranchida et al., 2018).

A rotationally invariant multi-scale formulation writes the total Hamiltonian as

$P(\sigma)\propto e^{-H/T}$3

where $P(\sigma)\propto e^{-H/T}$4 is harmonic, $P(\sigma)\propto e^{-H/T}$5 is Heisenberg exchange, $P(\sigma)\propto e^{-H/T}$6 is on-site uniaxial anisotropy defined from local neighbor geometry, and $P(\sigma)\propto e^{-H/T}$7 contains spin-lattice coupling terms such as the linear-in-displacement correction to exchange. Because each term depends only on scalar products among spins or on differences of position vectors, the Hamiltonian is translationally and rotationally invariant, and the continuous-time flow conserves total energy, total linear momentum $P(\sigma)\propto e^{-H/T}$8, and total angular momentum $P(\sigma)\propto e^{-H/T}$9 (Weißenhofer et al., 2022).

A distinct but related formulation is used in TSPIN, where the classical Lagrangian is extended by explicit spin kinetic terms and Nosé-Hoover chain variables. The canonical variables become $\mathbf r_i$0, $\mathbf r_i$1, and $\mathbf r_i$2, and the extended Hamiltonian is symplectic in NVE, NVT, and NPT form (Huang et al., 15 Jun 2025). This places spin and lattice variables in a unified Hamiltonian framework when machine-learning potentials are used for the energy $\mathbf r_i$3.

2. Splitting, rotation, and thermostat constructions

The dominant deterministic construction is symmetric operator splitting. In coupled spin-lattice dynamics, the Liouville operator is decomposed into drift, force-kick, and spin-precession parts. Representative second-order factorizations are

$\mathbf r_i$4

or, equivalently,

$\mathbf r_i$5

Each subflow is analytically integrable: drift updates positions, kick updates momenta, and the spin substep rotates each spin about its local effective field (Tranchida et al., 2018, Weißenhofer et al., 2022, Chen et al., 12 Jun 2026).

For single-spin precession, norm-preserving updates are written as exact or approximate orthogonal rotations. In the symplectic spin-lattice algorithm of Tranchida et al., the global spin propagator is further split into single-spin maps so that each update uses the most recent orientations of neighboring spins, and the Omelyan single-spin step preserves $\mathbf r_i$6 to the displayed order (Tranchida et al., 2018). In the rotationally invariant formulation of Weißenhofer et al., the spin subflow is $\mathbf r_i$7 with $\mathbf r_i$8, an exact rotation about the local field axis (Weißenhofer et al., 2022). In NEPSPIN, the spin step uses a midpoint iteration:

1. evaluate $\mathbf r_i$9,
2. predict $\mathbf p_i$0,
3. reevaluate the field,
4. average to a midpoint field $\mathbf p_i$1,
5. update $\mathbf p_i$2 (Chen et al., 12 Jun 2026).

A second construction replaces explicit lattice simulation by a conservative canonical thermostat. The Discrete Laplacian Thermostat (DLT) modifies spin dynamics according to

$\mathbf p_i$3

with the graph Laplacian

$\mathbf p_i$4

and noise built from the incidence matrix $\mathbf p_i$5 via $\mathbf p_i$6, where

$\mathbf p_i$7

Because $\mathbf p_i$8, the site noise covariance is $\mathbf p_i$9, which is the covariance required by the fluctuation-dissipation theorem to drive the system to $\mathbf s_i$0 (Cavagna et al., 2022).

A third construction targets the Landau-Lifshitz equation directly. The first-order method of 2026 combines a Gauss-Seidel predictor, a double-diffusion damping corrector, and a Crank-Nicolson preserving projection. In the undamped case, the last stage is

$\mathbf s_i$1

so the update matrix is orthogonal because $\mathbf s_i$2, implying $\mathbf s_i$3 (Xie et al., 11 Feb 2026).

At the spin-only end of the spectrum, the discrete-space-time anisotropic Landau-Lifshitz model gives an explicit discrete symplectic integration scheme in which the full step factorizes into odd and even two-body updates in a brick-wall circuit. The two-body map is a Poisson automorphism derived from a discrete zero-curvature representation, and the transfer matrix yields an infinite family of commuting conserved quantities (Krajnik et al., 2021).

3. Conservation laws, invariant measures, and long-time behavior

The defining feature of these methods is the exact or controlled preservation of specific structures at the discrete level. Exact spin-norm conservation is central. In spin-rotation substeps, each update is orthogonal, so $\mathbf s_i$4 or $\mathbf s_i$5 is preserved exactly or to machine precision (Tranchida et al., 2018, Chen et al., 12 Jun 2026). In the Landau-Lifshitz method, the Crank-Nicolson projection yields $\mathbf s_i$6 preservation by construction, and the reported maximum pointwise defect remains approximately $\mathbf s_i$7 (Xie et al., 11 Feb 2026).

Global conservation laws can be stronger. In DLT, the total magnetization $\mathbf s_i$8 is preserved exactly for every realization of the noise because the reversible term, the Laplacian drift, and the incidence-generated noise each sum to zero: $\mathbf s_i$9 The same construction yields exact sampling of the canonical ensemble because the Fokker-Planck operator has the unique stationary solution $\mathbf S_i$0 (Cavagna et al., 2022).

In rotationally invariant spin-lattice coupling, the split subflows $\mathbf S_i$1, $\mathbf S_i$2, and $\mathbf S_i$3 each preserve total linear momentum and total angular momentum, so the full symmetric Suzuki-Trotter composition preserves $\mathbf S_i$4 and $\mathbf S_i$5 exactly step by step. This is a discrete conservation result, not merely a continuous-time one (Weißenhofer et al., 2022).

Symplecticity alters the interpretation of energy conservation. In the coupled spin-lattice schemes of Tranchida et al., each submap is symplectic, and the symmetric product remains symplectic. The global error per step is $\mathbf S_i$6, so over a fixed interval the energy error grows like $\mathbf S_i$7, and long-time drift of first integrals remains bounded (Tranchida et al., 2018). TSPIN states the same principle through backward-error analysis: the symplectic and time-reversible map preserves a nearby modified Hamiltonian $\mathbf S_i$8, while the error in the actual extended Hamiltonian remains bounded by $\mathbf S_i$9 over exponentially long times (Huang et al., 15 Jun 2025). NEPSPIN likewise reports no secular drift in a pure Hamiltonian run and monitors a modified Hamiltonian $$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$0 with no runaway drift over millions of steps (Chen et al., 12 Jun 2026).

These distinctions matter because exact energy conservation is not universal. In DLT, exact energy conservation holds in the reversible limit $$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$1, where the method reduces to pure microcanonical spin dynamics; at finite $$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$2, the scheme becomes a conservative canonical dynamics that relaxes magnetic energy while preserving magnetization (Cavagna et al., 2022).

4. Representative algorithmic families

The recent literature contains several distinct families of structure-preserving integrators. They differ in whether the lattice is explicit or marginalized, whether the target ensemble is NVE or canonical, and which invariant is enforced exactly.

Method	Core mechanism	Preserved structure
Symplectic coupled SLD (Tranchida et al., 2018)	Suzuki-Trotter + single-spin updates	Spin norm; symplecticity
Rotationally invariant SLD (Weißenhofer et al., 2022)	$$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$3	Total $$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$4 and $$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$5
DLT (Cavagna et al., 2022)	Laplacian drift + incidence noise	Global magnetization; Gibbs measure
TSPIN (Huang et al., 15 Jun 2025)	Nosé-Hoover-chain Hamiltonian + Strang splitting	Symplectic NVE/NVT/NPT dynamics
NEPSPIN (Chen et al., 12 Jun 2026)	Verlet-type $$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$6-$$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$7 splitting + exact spin rotations	Spin norm; time-reversibility
Landau-Lifshitz method (Xie et al., 11 Feb 2026)	Gauss-Seidel + double diffusion + Crank-Nicolson	Norm preservation; unconditional stability in implicit form

The coupled spin-lattice symplectic algorithm implemented in LAMMPS is explicitly designed for large spin-lattice systems and combines Suzuki-Trotter decomposition for non-commuting variables with a sectoring strategy for parallel domain decomposition (Tranchida et al., 2018). The rotationally invariant formalism corrects earlier Suzuki-Trotter decompositions and embeds exchange, anisotropy, and spin-lattice coupling in a Hamiltonian that is explicitly translationally and rotationally invariant (Weißenhofer et al., 2022). DLT occupies a different niche: it is a low-cost, structure-preserving canonical integrator for spin systems with a conserved order parameter and does not actually simulate the lattice, even though its parameter $$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$8 is quantitatively connected to microscopic spin-lattice couplings (Cavagna et al., 2022).

TSPIN and NEPSPIN extend the same structure-preserving logic to machine-learning potentials. TSPIN introduces explicit spin kinetic terms and thermostat variables in a unified Hamiltonian, yielding second-order, time-reversible, symplectic splitting for NVE, NVT, and NPT ensembles (Huang et al., 15 Jun 2025). NEPSPIN combines a structure-preserving spin-lattice integrator with a spin-constrained density-functional-theory-trained neuro-evolution potential and augments the algorithm with fused force-torque kernels, SVE2 vectorization, and SME outer-product acceleration (Chen et al., 12 Jun 2026).

5. Numerical validation and application domains

The most detailed validation of a conservative thermostat appears in DLT on the 3D Heisenberg antiferromagnet with

$$\mathcal H(\{\mathbf r_i,\mathbf p_i,\mathbf S_i\}) = \sum_{i=1}^N \frac{|\mathbf p_i|^2}{2m_i} + E\bigl(\{\mathbf r_i\},\{\mathbf S_i\}\bigr),$$9

on an $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$0 periodic cubic lattice, $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$1, with $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$2. Starting from $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$3 initial states, the spin energy $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$4 relaxes exponentially to its canonical value on a time scale $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$5. The staggered magnetization $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$6 vanishes at $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$7. Finite-size scaling of the susceptibility at its peak gives $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$8 with fitted $$\dot{\mathbf r}_i=\frac{\mathbf p_i}{m_i},\qquad \dot{\mathbf p}_i=-\nabla_{r_i}E,\qquad \dot{\mathbf S}_i=\mathbf S_i\times \mathbf b_i,\qquad \mathbf b_i=-\nabla_{S_i}E.$$9, compared with the exact 3D Heisenberg value $\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$0. At criticality, the lowest-mode relaxation time scales as $\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$1 with $\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$2, in precise agreement with the exact Model G result $\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$3. For $\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$4, the transverse dynamical structure factor shows two symmetric peaks whose dispersion matches

$\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$5

with no fitting parameters. Increasing $\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$6 broadens and shifts the spin-wave peaks downward, reproducing magnon-phonon damping observed in SD+MD; typical estimates in magnetic insulators give $\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$7 in the range $\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$8–$\mathcal H_{\rm mag}(\mathbf r,\mathbf s)$9 (Cavagna et al., 2022).

For explicitly coupled dynamics, the parallel symplectic algorithm was tested on 500–2000 cobalt atoms in fcc geometry, where both total energy and the norm of the total magnetization fluctuate only within bounds proportional to $V(r_{ij})$0, confirming second-order accuracy (Tranchida et al., 2018). In the rotationally invariant multi-scale framework, simulations of a ferromagnetic nanoparticle recover not only the ferromagnetic resonance mode but also another low-frequency mechanical response and a rotation of the particle according to the Einstein-de-Haas effect (Weißenhofer et al., 2022).

Machine-learning-based integrators extend validation to modern atomistic workloads. For FCC Fe at $V(r_{ij})$1 and $V(r_{ij})$2, TSPIN in NVT and NPT ensembles reports $V(r_{ij})$3 for time steps $V(r_{ij})$4, $V(r_{ij})$5, and $V(r_{ij})$6 over $V(r_{ij})$7, whereas the conventional LLG method exhibits drifts one-two orders of magnitude larger. On an NVIDIA V100 GPU, classical MD scales as linear $V(r_{ij})$8, TSPIN SLD is also $V(r_{ij})$9 and nearly indistinguishable from MD, while LLG-SLD and MD+MC are $P(\sigma)\propto e^{-H/T}$00 (Huang et al., 15 Jun 2025).

NEPSPIN pushes the same structure-preserving philosophy to extreme scale. Deployed on the LineShine exascale supercomputer, the full application scales to $P(\sigma)\propto e^{-H/T}$01 million CPU cores with $P(\sigma)\propto e^{-H/T}$02 weak-scaling efficiency, enabling simulations of $P(\sigma)\propto e^{-H/T}$03 trillion atoms and an equal number of spins while reaching $P(\sigma)\propto e^{-H/T}$04 PFLOPS in double precision. The paper reports a seven orders-of-magnitude speedup over prior spin-aware methods and states that the simulations directly resolve real-temperature skyrmion nucleation and reorganization at previously inaccessible scales (Chen et al., 12 Jun 2026).

The Landau-Lifshitz structure-preserving method illustrates a different validation regime. In one-dimensional tests with exact solution and $P(\sigma)\propto e^{-H/T}$05, it shows first-order convergence in time, second-order convergence in space, and norm preservation to machine precision, with $P(\sigma)\propto e^{-H/T}$06. In three-dimensional tests with $P(\sigma)\propto e^{-H/T}$07, the same first-order-in-time and second-order-in-space behavior is reported. The fully implicit versions require no upper bound on $P(\sigma)\propto e^{-H/T}$08 for stability and permit $P(\sigma)\propto e^{-H/T}$09–$P(\sigma)\propto e^{-H/T}$10 orders of magnitude larger $P(\sigma)\propto e^{-H/T}$11 than a standard explicit solver (Xie et al., 11 Feb 2026).

6. Conceptual scope and common distinctions

A recurring distinction is between explicit spin-lattice dynamics and effective spin-only dynamics informed by lattice physics. DLT belongs to the latter category: it turns microcanonical spin dynamics into a conservative canonical dynamics, preserves the constants of motion, and relaxes magnetic energy without actually simulating the lattice. Its additional parameter obeys

$P(\sigma)\propto e^{-H/T}$12

in the Markovian limit of an explicit spin-phonon model, where $P(\sigma)\propto e^{-H/T}$13 is the spin-lattice coupling and $P(\sigma)\propto e^{-H/T}$14 is the lattice stiffness (Cavagna et al., 2022). This makes $P(\sigma)\propto e^{-H/T}$15 a quantitatively connected surrogate for microscopic spin-lattice coupling rather than an arbitrary damping constant.

Another common distinction is that structure preservation does not always mean the same invariant. In symplectic Hamiltonian splitting, the characteristic result is bounded energy error and preservation of a modified Hamiltonian, not exact total energy at finite $P(\sigma)\propto e^{-H/T}$16 (Tranchida et al., 2018, Huang et al., 15 Jun 2025, Chen et al., 12 Jun 2026). In rotationally invariant splitting, the central exact discrete invariants are total linear and angular momentum (Weißenhofer et al., 2022). In conservative stochastic thermostats, the exact target is the canonical stationary distribution together with conserved magnetization (Cavagna et al., 2022). In norm-preserving Landau-Lifshitz solvers, the primary exact discrete invariant is $P(\sigma)\propto e^{-H/T}$17, while energy-stability is argued through dissipativity of the diffusion substep and energy-conserving character of the precessional Crank-Nicolson stage (Xie et al., 11 Feb 2026).

The literature also distinguishes between general-purpose geometric integrators and integrable discrete-time spin maps. The anisotropic lattice Landau-Lifshitz model in discrete space-time is an explicit discrete symplectic integration scheme built from a discrete zero-curvature representation; it preserves the Poisson structure exactly and generates commuting local integrals from the transfer matrix (Krajnik et al., 2021). This suggests a broader methodological continuum: from exact integrable spin discretizations, through symplectic and momentum-preserving spin-lattice splittings, to conservative thermostats and machine-learning-based Hamiltonian formulations.

A further implication drawn explicitly in the machine-learning setting is that non-symplectic integration can produce poor energy conservation and excessive computational costs. TSPIN addresses this by treating spins and lattice simultaneously within a symplectic Hamiltonian framework, while NEPSPIN couples a structure-preserving integrator to a machine-learned spin-lattice potential and architecture-specific optimization. Taken together, these developments show that the phrase “structure-preserving spin-lattice integrator” refers less to a single update formula than to a design principle: the numerical map is built so that the invariants, symmetries, and ensemble structure of the underlying spin-lattice model survive discretization as exactly as the chosen formulation permits (Huang et al., 15 Jun 2025, Chen et al., 12 Jun 2026).