Anisotropic Landau-Lifshitz Model in Discrete Space-Time

We construct an integrable lattice model of classical interacting spins in discrete space-time, representing a discrete-time analogue of the lattice Landau-Lifshitz ferromagnet with uniaxial anisotropy. As an application we use this explicit discrete symplectic integration scheme to compute the spin Drude weight and diffusion constant as functions of anisotropy and chemical potential. We demonstrate qualitatively different behavior in the easy-axis and the easy-plane regimes in the non-magnetized sector. Upon approaching the isotropic point we also find an algebraic divergence of the diffusion constant, signaling a crossover to spin superdiffusion.


Introduction
The advent of powerful computational tools [1,2] has tremendously advanced our understanding of many-body physics in and out of equilibrium [3][4][5][6]. Dealing with thermodynamic systems of strongly interacting degrees of freedom nonetheless still presents a dauting task in general. In order to efficiently simulate many-body quantum dynamics on a computer one typically confronts the issue of growing complexity which arises when distant regions of space become entangled. While classical many-body dynamical systems on the other hand do not suffer this issue, computing time-evolution of statistical ensembles often proves to be similarly challenging and require an immense amount of computational resources.
Amidst recent experimental breakthroughs with cold-atom technologies [7,8], there has been an upsurge of interest in the study of out-of-equilibrium phenomena. Nowadays, optical lattices not only enable manufacturing of tailor-made nanodevices for processing quantum information, but also offer a versatile platform to directly probe the relaxation process at finite energy density in the systems both near and far from equilibrium [9][10][11][12].
One major point of worry concerning numerical simulations of integrable dynamics is that integration schemes based on naïve (e.g. Trotter-Suzuki, or Runge-Kutta) discretizations invariably destroy integrability. On sufficiently short time scales this should not be a major concern. By contrast, even weakly broken integrability is likely to induce spurious effects at late times and thus preclude reliable extraction of transport coefficients or dynamical exponents. In addition statistical field theories are commonly plagued by UV divergences (even when the local target space is a compact manifold). These drawbacks can both be obviated in a fully discrete setting.
Integrable many-body dynamical systems in discrete time [53,54] and classical cellular automata [43,55] have recently attracted much interest. An integrable space-time discretization of the isotropic classical Heisenberg ferromagnet has been recently obtained in Ref. [56], and subsequently generalized to a large class of 'matrix models' [54] which are globally invariant under the action of simple Lie groups. In this work, we construct an anisotropic deformation of the classical SU (2) ferromagnet obtained in Ref. [56], representing a 'brick-wall type' circuit composed from elementary two-body symplectic maps. Our model can be thus regarded as an 'integrable Trotterization' [54,[56][57][58] of the anisotropic lattice Landau-Lifshitz field theory [59][60][61] (the classical counterpart of the celebrated quantum XXZ spin chain). A fully discrete integrable Landau-Lifshitz equation has, to the best of our knowledge, not been systems (see e.g. Refs. [61,65,66]). To realize parallel transport on the auxiliary lightcone lattice we introduce the lattice shift operators [67] along the light-cone directions (i.e. characteristics ± t = const) φ n+1,m = L (+) n,m (λ)φ n,m , φ n,m+1 = L (−) n,m (λ)φ n,m . (2.2) Here the local 'matrix propagators' L (±) represent certain matrix-valued functionals of physical variables which, in addition, depend analytically on a free complex parameter λ. In order to ensure consistency of such parallel transport, there should be no ambiguity in the order of propagation when passing from φ n,m to φ n+1,m+1 . In this case, the shift matrices L (±) constitute the Lax pair and λ is called the spectral parameter. More generally, one can also allow the spectral parameter to have additional dependence on local 'inhomogeneities' along an initial sawtooth. We shall make a special 'homogeneous choice' by requesting L ± ≡ L(λ ± ; S), i.e. using a pair of spectral parameters λ ± depending only on the light-cone direction (but not on position or time t).
Commutativity of the light-cone shifts is neatly encapsulated by the following discrete zerocurvature property [68,69] around an elementary square plaquette of the light-cone lattice, 1 L(λ + ; S 2 )L(λ − ; S 1 ) = L(λ − ; S 2 )L(λ + ; S 1 ). (2.3) where we have assumed, for notational convenience, that the 'input' variables S 1,2 sit at lattice sites 1 and 2. Notice that each of the light-cone Lax matrices L(λ ± ) depend on the local 'edge variable' S on the corresponding edge, whereas primed variables S are being used as a shorthand notation for the propagated spins, i.e. original variables S shifted by one unit in the time direction (as depicted in Figure 1). The discrete curvature is fulfiled if and only if the time-updated 'output' variables S 1,2 are appropriately linked to 'input' variables S 1,2 . This requirement allows us to interpret the flatness condition as an (implicit) specification of a local time-propagator, i.e. a local two-body map Φ : (S 1 , S 2 ) → (S 1 , S 2 ) for any spatially adjacent pair of sites.
The zero-curvature condition implies the integrability of the model through a property called isospectrality. As explained in Appendix C, isospectrality gives rise to an infinite set of local conserved quantities in involution.
Many-body propagator. Let us for the time being leave the two-body propagator Φ defined implicitly via the zero-curvature property. In the following, we first outline how to use Φ as the elementary building block of a many-body dynamical map in the form of a 'brick-wall circuit'. To this end, let M 1 ∼ = S 2 denote a local phase space (a 2-sphere) of a canonical spin. The full phase space M L of an L-site chain is thus given simply by the L-fold Cartesian product, M L ≡ M ×L 1 , and accordingly we introduce a many-body dynamical map Φ full : M L → M L .
By virtue of the light-cone geometry, the full time-dynamics consists of an alterating sequence of 'odd' and 'even' local propagators, respectively. This results in a staggered circuit as depicted in Figure 2. By making use of the embedding prescription, for j = 1, . . . , L − 1, where I : M 1 → M 1 stands for a local (single-site) unit map I(S) ≡ S, (Φ (L) correspondingly affects the Lth and the first spin), the full propagator Φ full for two units of time t → t + 2 decomposes as with odd and even propagators further factorizing into commuting two-body maps, So far our construction has been entirely formal. In what follows, we derive the map Φ full that corresponds to an integrable space-time discretization of the anisotropic lattice Landau-Lifshitz model, as introduced by Sklyanin [60,61,70]. In this view, Φ full represents its 'integrable Trotterrization'.

Sklyanin Lax matrix
In this section, we construct an integrable time-discrete analogue of the anisotropic lattice Landau-Lifshitz model. To facilitate the computations, instead of canonical spins S, it proves more convenient to employ an axially 'deformed spin' with components {K, S + , S − }, enclosing a Poisson algebra This Poisson algebra admits a quadric (Casimir) function satisfying {C 0 , S a } = 0. In terms of 'Sklyanin spins' 2 , with components {S a } 3 a=0 , with the Hamiltonian of the lattice Landau-Lifshitz model takes a compact form with couplings g 0 = sinh 2 ( ), g 1 = g 2 = 1 and g 3 = cosh 2 ( ). We note that A q is just the trigonometric limit of 'elliptic spins' that constitute the socalled quadratic Poisson-Sklyanin algebra, cf. Appendix A for further information. The algebra A q becomes non-degenerate on symplectic leaves on which C 0 takes a constant value c 0 . We consider two physically distinct regimes: (i) the easy-axis regime, with a real anisotropy (or interaction) parameter ∈ R + , 2 Classical canonical spins can be recovered in the → 0 limit by simultaneously rescaling S a → S a / for a ∈ {1, 2, 3}.
Our construction, which we outline in turn, applies to both of these regimes. Before that, we first briefly discuss the symplectic structure of the local space space. To obtain the easy-axis lattice Landau-Lifshitz model, the Casimir function C 0 has to be fixed to c 0 = sinh 2 ( ). (2.12) In this way, we select a two-dimensional non-degenerate Poisson submanifold M ea inside R 3 which is diffeomorphic to a 2-sphere. In the easy-plane regime, obtained under analytic continuation → iγ, the submanifold M ep is also topologically equivalent to a 2-sphere provided γ ∈ [−π/2, π/2]. The upshot here is that there exist a smooth bijective mapping to the symplectic sphere S 2 , implying that variables K and S ± can be parameterized in terms of classical canonical spins S. A local mapping from S 2 to M ea is given by [61] K = e S 3 , the form of which follows from the Casimir function, see Eqs. (2.9),(2.12). Two imporatnt remarks are in order at this stage: (I) in the easy-plane regime (with → iγ and F iγ (s)), Eq. (2.13) provides a bijective map only inside a compact interval of anisotropies |γ| ≤ π/2; (II) there exist other Poisson submanifolds (classified in [61,73]) where the classical Sklyanin algebra becomes non-degenerate. Those are not expressible in terms of canonical spins and thus we do not consider them here.
The discrete zero-curvature condition (2.3) may be formally viewed as a matrix refactorization problem. The task ahead of us is to find its explicit solution, i.e. expressing the output (primed) variables in terms of the input spin variables. We achieve this by first casting the Sklyanin Lax matrix in the form where we have for convenience introduced a multiplicative spectral parameter z ∈ C, In Appendix A we detail out how the above Sklyanin Lax matrix arises from the limit of the quantum Lax operator associated to the quantum Yang-Baxter algebra. The task of refactorization is unfortunately not straightforward. To begin with, Eq. (2.3) merely provides an implicit rule for the time-propagated input variables, and there is no obvious way of recasting it as an explicit map Φ. 4 This difficulty can be elegantly overcome by exploting certain factorization properties of the Lax matrix into more fundamental constituents.

Factorization
In this section, we outline how the Lax matrix can be further decomposed into more elementary constituents. Such factorization properties are indeed deepely rooted in the algebraic properties of certain types of Hopf algebras that are associated to simple Lie algebras or deformations thereof. For instance, it is well known that the universal quantum R-matrix admits a triangular (Borel) decomposition [74][75][76][77]. The resulting 'partonic' Lax operators can be used to construct the Baxter Q-operators [78]. There exist various equivalent constructions in the literature, see e.g. Refs. [79][80][81] and Refs. [82][83][84]. We in turn demonstrate how the entire construction quite naturally elevates to the classical level by considering its semiclassical limit (outlined in Appendix A). Below we collect the relevant formulae without delving too much into formal considerations.

Weyl variables.
To facilitate the factorization of the Lax matrix, we parametrize Sklyanin variables in terms of a classical Weyl-Poisson algebra. The latter is spanned by a pair of generators, x, y, satisfying the multiplicative canonical Poisson bracket {x, y} = i x y, (2.17) via the prescription where we have simultaneously defined, for compactness of notation, ν ≡ e . We accordingly make the identification y = ν −1 K = e (S 3 −1) . (2.19) This implies that y is real and positive in the easy-axis case, whereas in the easy-plane regime it becomes unimodular. In Weyl variables, however, the reality condition S + = S − is not manifestly satisfied. To enforce it, the variable x has to be restricted to a 'physical' submanifold by demanding the squared modulus to obey (2.20)

Factorization. By introducing
along with a pair of spectral parameters, satisfying uv = ν 2 , u/v = z 2 , the Lax operator from Eq. (2.15) can be factorized as in terms of diagonal 'coordinate' matrices X and Y 24) and the 'spectral' matrices U and V Elementary exchange relations. We next describe various elementary exchange procedures involving L x,y (U, V ). Firstly, interchanging the spectral matrices specifies a map (x, y) → (x , y ). The canonical Weyl relations are preserved 6 by additionally assuming y = y, in which case Secondly, there exist quadratic exchange relations of the type which are canonical provided that From Eq. (2.27) we can deduce the following transformation of (x, y) → (x , y = y):

33)
∈ {1, 2}, while analogously the exchange of U 1 and U 2 according to Eq. (2.31) yields the same expression with v replaced by u : By the same reasoning as previously (only swapping the roles of x and y ), transformations (2.33) and (2.34) are both canonical.
Composite exchange relations. By using the above elementary exchange relations it is possible to also exchange U 1 with U 2 between a pair of Lax matrices This can be achieved by a composition of three elementary exchanges In an analogous fashion we can exchange via the following sequence of elementary exchanges Solving the zero-curvature condition. To obtain the two-body propagator Φ the discrete zero-curvature condition, Eq. (2.3), Their composition gives an explicit propagator Φ that solves the zero-curvature condition (2.39). Because the elementary exchanges are canonical transformations, the resulting propagator is automatically canonical. Explicit formulas for these compositions are rather cumbersome and we display them in Appendix B.

Two-body propagator
In analogy to the isotropic case [56], the discrete 'Trotter' time-step τ ∈ R enters the construction through the difference of the additive alternating (i.e. staggered) spectral parameters λ ± ≡ λ ± τ /2, carried by the light-cone Lax matrices, cf. Eq. (2.3). In terms of spectral parameter u and v, this amounts to setting recalling that ν = e . With the above prescription, the exchange relations (given by equations (2.35) and (2.37)) and the propagator Φ simplify considerably. Explicit expressions can be found in Appendix B).
The elementary propagator Φ τ, can alternatively be written explicitly using Sklyanin variables K, S ± . With aid of Eqs. (2.18) and some exercise, we find (where w = e i τ ) with We remind the reader that in the easy-axis regime K ∈ R + , whereas in the easy-plane regime γ S 3 = arg(K) ∈ [−π/2, π/2]. Therefore, there is no ambiguity in Eqs. (2.42) when taking the square roots to determine K .

Properties.
There are a number of important properties worth highlighting: • In the continuous-time limt τ → 0, the propagator Φ τ, reduces to the identity map.
• The mapping (2.42) manifestly preserves the product K 1 K 2 , implying that the local propagator conserves the third component of the spin, This local U (1) symmetry directly implies conservation of S 3 tot ≡ L =1 S 3 under the full propagator Φ full which, in the Hamiltonian limit, τ → 0, becomes the Noether charge of the uniaxially anisotropic Landau-Lifshitz model.
• In the easy-axis regime, the local propagator (2.42) is a periodic function of interaction , with period 2π/τ .
• The two-body propagator Φ τ, is invariant under the transformation It is therefore sufficient to consider only the interval ∈ [0, π/τ ].
• In the easy-plane regime, in the limit τ → ∞ with γ fixed (i.e. |w| → ∞), Φ τ,iγ reduces to a permutation lim • The two-body propagator Φ is a classical Yang-Baxter map [67,[85][86][87], obeying the braided Yang-Baxter relation on the Cartesian product M ×3 1 (suppressing the anisotropy parameter) Φ (1,2) where Φ (i,j) signifies that Φ operates non-identically on i-th and j-th copy of M ×3 1 . Figure 3: Discrete Noether current, shown at an even site 2 : the time-increment t → t + 2 of charge density, q t+2 2 − q t 2 /τ equals the difference of two nearby current densities at consecutive times, j 2t+1 2 +1 −j 2t 2 . Current densities j t , depicted by pink arrows, can be compactly expressed as the difference of charge densities at consecutive times, j t = (q t+1 − q t )/τ . Discrete Noether current. The lattice Landau-Lifshitz model, cf. Eq. (2.11), is invariant under global rotations about the third axis. This imples, by virtue of the Noether theorem, a globally conserved Noether charge (magnetization) in the system. We have just explained, cf. Eq. (2.44), that the same property holds even in the discrete-time setting.
By virtue of the global U (1) symmetry, the dynamics admits a conserved Noether current. Its spatial component, i.e. the 'charge density', is given by the projection of S onto the third axis, q = S 3 . The associated temporal component is the 'current density' which we subsequently denote by j. Owing to the even-odd staggered structure of the full propagator Φ full , the discrete continuity equations at even and odd sites are of the form respectively. While the charge densities q t ≡ q( , t) are ultralocal, namely they sit at site and time t, the associated current density corresponds to their forward differences 7 which can be immediately verified by plugging the expression into the discrete continuity equations (2.49) and using the conservation property (2.44). Notice that q t+1 is a function of two adjacent charge densities at the previous time slice, namely q t and q t ±1 , cf. Figure 3, that enter the local propagator Φ.
Isotropic, continuous time and field-theory limits. The two-parametric dynamical map Φ full τ, admits various important limits: (i) the limit of vanishing anisotropy → 0, (ii) the continuous-time limit τ → 0 and (iii) the continuous space-time (i.e. field-theory) limit by additionally taking the long wavelength limit. The resulting models, outlined in Appendix D, can be all regarded as members of the 'Landau-Lifshitz family'.

Magnetization transport
We demonstrate the utility of the constructed dynamical system by studying magnetization transport in thermal equilibrium. We focus on numerical computation of the linear transport coefficients that quantify magnetization transport on a large spatio-temporal (i.e. hydrodynamic) scale, namely the spin Drude weight and the spin diffusion constant. After introducing the equilibrium measure and transport coefficients we proceed by systematically analyzing the easy-axis and the easy-plane regimes.

Transport coefficients
Equilibrium ensemble. We are interested in the dynamical response in our model at finite magnetization density. To this end, we introduce a one-parametric family of local measures which can be interpreted as the grand-canonical measure with chemical potential µ. 8 We accordingly define a separable stationary measure with chemical potential µ enforcing a non-zero average of magnetization. Invariance of ρ L under the full dynamics Φ full follows from the fact that the quadratic Sklyanin bracket is preserved by the action of Φ (see Ref. [54] for an analogous formal proof). We now pass over to the thermodynamic limit by sending L → ∞. The grand-canonical free energy per site is given by where represents the onsite grand-canonical partition function, with dΩ denoting the volume element on S 2 . Writting dΩ full ≡ L =1 dΩ(S ), the corresponding thermodynamic average of any local observable O is therefore computed as For example, the average of magnetization density reads In the following we accordingly take L → ∞ and operate in the thermodynamic setting.

Conductivity.
In the linear-reponse regime, transport coefficients are encoded in the real part of the frequency-dependent conductivity σ(ω), which is in general decomposed as The δ-peak owes its existence to the presence of long-lived excitations which ballistically transport the charges through the system, characterized by a coefficient D called the Drude weight. The d.c. component of the 'regular' part, that is σ reg (0), characterizes transport on sub-ballistic scales. Provided σ reg (0) is finite, it can be associated to the spin diffusion constant defined as where χ denotes static spin susceptibility, In integrable interacting systems, the diffusion constant quantifies the diffusive spreading which quasiparticle excitations exhibit upon ellastic collisions (see e.g. Refs. [41,88]). As customary, we define the dynamical structure factor with a sum rule ∈Z S( , t) = χ. The spin Drude weight governs the time-asymptotic growth of its second moment (3.10) where the subleading correction, which equals 2χ D, stores information about the spin diffusion constant D. We note that numerically extracting D from the O(t) broadening of the second moment (3.10) has not proven particularly reliable.
We shall employ an alternative approach and define the (linear) transport coefficients through the Kubo formula. The extensive magnetization current J at physical times 2t follows directly from the discrete continuity equation (2.49) and is given by The spin Drude weight corresponds to the asymptotic value of the current autocorrelation function whereas the spin diffusion constant is the integrated current autocorrelator with the asymptotic value subtracted Lastly, we introduce the dynamical exponent α that characterizes the late-time decay of the dynamical spin structure factor S(0, t) ∼ t −α . (3.14) Ballistic spreading corresponds to exponent α = 1. On the other hand, we can speak of normal diffusion when (i) α = 1/2 and (ii) the dynamical structure factor converges asymptotically to a Gaussian stationary scaling profile,

Numerical simulations: easy-axis regime
We first focus our analysis on the spin Drude weight and diffusion in the easy-axis regime.
The results are shown in Figure 4. The asymptotic stationary profiles and the numerically estimated dynamical exponents are shown in Figure 5.
Drude weight. The spin Drude weight is finite for all anisotropies and arbitrary finite magnetization density µ > 0. In the unmagnetized ensemble (µ = 0), corresponding to a uniform measure ρ 1 = 1/Vol(S 2 ) = 1/4π, the Drude weight vanishes. This point is distinguished by the discrete Z 2 symmetry, corresponding to the global reversal of all spins, S ± → S ∓ , S z → −S z , which is a canonical transformation, i.e. preserving the Poisson bracket. We find, moreover, that both transport coefficients vanish when → π/τ , consistent with the dynamics becoming trivial in that limit, cf. Eq. (2.46). Fixing anisotropy to = 1, we find a non-monotonic dependence of D on the chemical potential µ, see the inset plot in Figure 4, which is consistent with the picture of dressed quasiparticles [41,42]: in an unpolarized ensemble, the dressed magnetization vanishes and so does the spin Drude weight. For small µ, the quasiparticles behave paramagnetically, with dressed magnetization behaving as ∼ µ. For strong polarizations, µ → ∞, one approaches a Z 2 -degenerate quasiparticle pseudovacuum where the quasiparticles become effectively free.
In this regime, however, the quasiparticle density vanishes and the spin Drude weight drops to zero.

Diffusion constant.
The spin diffusion constant D is finite for all µ > 0. We find that it decreases monotonically with increasing µ. At any given anisotropy , D grows upon lowering µ, and attains its maximal value at µ = 0. There, upon approaching the isotropic point, → 0, the diffusion constant D diverges in (approximately) algebraic fashion, D ∼ 1/ κea , as shown in Figure 4 (right inset). This conforms with a theoretical expectation κ ea = 1 based on the known behavior in the quantum Heisenberg chain [35]. 9 A reliable estimation of κ ea from the numerical data, as shown in Figure 4 (b, inset, dashed black line), is however rather difficult due to a divergent relaxation timescale in the proximity of the isotropic point. Singularity of the spin diffusion constant signals the onset of superdiffusion [89] which has recently been under intense scrutiny [31,33,54,56,72,[90][91][92]. There exists ample numerical evidence [33,34,54,56,72,93,94] that this phenomenon belongs to the universality class of the Kardar-Parisi-Zhang (KPZ) equation [95]. In Figure 5, we numerically corroborate the predicted anomalous dynamical exponent exponent α = 2/3. While the dynamics is ballistic for all values of > 0 and µ > 0, the finite-time dynamical exponent estimated based on (3.14) is found to be ballistic, α = 1, only below a threshold value of (which increases with increasing µ). This is due to the fact that the value of D is small in comparison with D. By contrast, at µ = 0, the estimated exponent is approximately diffusive, α ≈ 1/2 for 1, while when approaching → 0 it starts to drift towards the expected superdiffusive exponent α = 2/3. At µ = 0, the stationary cross-sections shown in Figure 5

Numerical simulations: easy-plane regime
We now repeat the above numerical analysis for the easy-plane regime. The spin Drude weights and diffusion constants are displayed in Figure 7, while the asymptotic stationary profiles of the dynamical structure factors, alongside the finite-time dynamical exponents, are shown in Figure 8 and Figure 9.
Drude weight. In the easy-plane regime, the spin Drude weight D takes a non-zero value for all values of µ and γ, except at the isotropic point at zero magnetization (µ = γ = 0). Upon lowering γ, the spin Drude weight D now monotonously decreases, for any fixed µ, to a limiting value at γ = 0 which is non-zero for µ = 0. It is worth highlighting that D does not vanish even in a non-magnetized ensemble, in stark contrast with the easy-axis regime. Altough this is perfectly aligned with previous numerical studies of the lattice Landau-Lifshitz model [71,72], a proper theoretical explanation is currently still lacking. At this juncture, we can offer two complementary interpretations: (a) from the quasiparticle point of view, finite spin Drude weight in a non-magnetized sector indicates that certain excitations attain finite dressed magnetization even in the µ → 0 limit, mirroring the situation in the easy-plane phase of the quantum Heisenberg XXZ spin chain [42,96], (b) from the viewpoint of the Mazur-Suzuki hydrodynamic projection [97,98] (for more information see e.g. [88,99,100]) a finite value of D necessitates the existence of local (or quasilocal, see Refs. [99,101,102]) conserved quantities with finite spin-charge susceptibility even at µ → 0, that is J Q c 0 > 0. Conserved quantities with a finite overlap with spin current have indeed been found in the gapless phase of the quantum Heisenberg chain (see Refs. [103][104][105][106] and [102] for a review). They are commonly referred to as the 'Z-charges'. Presently however, the standard local conservation laws Q ± k;τ extracted by expanding the trace of the monodromy matrix as an analytic series (see the derivation in Appendix C) all satisfy J Q ± k;τ c 0 = 0. This indeed follows directly from a simple symmetry argument: while all local charges Q k;τ remain invariant under the spin-reversal transformation, the spin current J flips a sign.
Diffusion constant. The diffusion constant D is finite for all values of γ provided µ > 0. In the µ → 0 limit however, it once again diverges as D ∼ 1/γ κep algebraically with an estimated exponent κ ep = 1, cf. Figure 7 (b, inset). We can observe, see Figure 8 (a, blue line), a crossover behavior in the vicinity of the isotropic point from ballistic to superdiffusive dynamics. In the easy-plane regime, the stationary profiles ( Fig. 9) have two peaks which are indicative of ballistic spreading. At large τ (with fixed γ) the two-particle map reduces to the permutation, see Eq. (2.47), meaning that the dynamics becomes non-interacting and hence the structure factor concentrates near the light-cone boundary peaked at /t = ±2.

Conclusion
We have introduced a new integrable model in discrete space-time, representing classical spins which interact locally via an anisotropic spin-spin interaction. The model can be perceived as an integrable discrete-time analogue of the uniaxailly anisotropic Landau-Lifshitz field theory. The key ingredient of our algebraic construction is the discrete version of the zerocurvature condition on an auxiliary light-cone lattice, which implicitly defines the elementary time-propagator. We managed to express the elementary propagator explicitly as a two-body symplectic transformation by exploiting certain factorization properties of the Lax matrix. The full propagator takes the form of a classical circuit made of symplectic maps. The two free parameters of the model are the Trotter time-step and the interaction anisotropy, the latter can be chosen real or purely imaginary (corresponding to the easy-axis and easy-plane regimes, respectively).
Our discrete model provides an efficient explicit integration scheme that manifestly preserves integrability. As an important immediate physics application, we numerically studied magnetization transport in grand-canonical stationary ensembles at different values of magnetization density. With only a moderate computational effort, we compute the linear transport coefficients (the spin Drude weight and spin diffusion constant) with aid of the Kubo formula, both in the easy-axis and easy-plane regimes.
The spin Drude weight is a non-negative quantity. As expected, we found that it only vanishes in (a) the limit of strong polarization and in (b) the easy-axis regime at zero net magnetization. On sub-ballistic scales, there is a finite diffusive correction characterized by the spin diffusion constant. According to expectations, we have found that the spin diffusion constant blows up when approaching the isotropic point (within either regime). This divergence is numerically estimated to be roughly inversely proportionally to the anisotropy strength. Close to the isotropic point, we encountered discernible finite-time crossover effects which signal the onset of spin superdiffusion [42] (characterized by the dynamical exponent α = 2/3).
In the easy-plane regime, the spin Drude weight retains a finite value even at zero net magnetization, consistent with the findings of Refs. [71,72]. It is well-known indeed that the same behavior occurs in the the gapless phase of the Heisenberg quantum spin chain, where the non-vanishing of the spin Drude weight is intimately linked with the existence of the so-called quasilocal conservation laws [102][103][104] of odd parity under the spin-reversal transformation. The most suggestive explanation of this peculiarity is that there exist additional 'hidden' (quasi)local conservation laws in the classical spin model as well -a classical analogue of the Z-charges.
The hope is that our numerical data can serve as an accurate test of these theoretical predictions. The spin Drude weight and diffusion constant are in principle amenable to analytic computation, e.g. using the universal formulae for the Drude weights [107,108] and diffusion constants [109][110][111] obtained within the formalism of generalized hydrodynamics [28,29]. This computation however requires the knowledge of the quasiparticle content of the model and the associated thermodynamic state functions [88], which are yet to be derived, e.g. by means of the inverse scattering techiques [61,65,66,112] adapted to the fully discrete setting. We leave this task for future research.

Appendices A Semiclassical limit of quantum algebra
The algebraic structures pertaining to the lattice Landau-Lifshitz model can be derived systematically by taking the semiclassical limit of a quasi-triangular Hopf algebra (i.e. the 'quantum group') U q (su 2 ). For the reader's convenience, we summarize the key steps of the derivation.

Fundamental commutation relation.
A natural starting point is to consideer the socalled fundamental commutation (or RLL) relation of the difference form [70,74,113] representing an operatorial identity on the tensor product of three vector spaces of the form V 1/2 ⊗ V 1/2 ⊗ V S . We have employed the standard embedding prescription in which the subscript indices specify on which spaces the operators act non-identically. Lax matrices act on C 2 ⊗ V S , where V S pertains to a unitary irreducible representation of 'Sklynanin quantum spins'Ŝ a (with a ∈ {0, 1, 2, 3}) of dimension 2S + 1. The Lax operator assumes the form Quantum R-matrix. The R-matrix acts on two copies of V 1/2 ∼ = C 2 which are regarded as auxiliary fundamental spins (subsequently represented by Pauli matrices σ a ). The R-matrix swaps (intertwines) the spectral parameters of the Lax operators.
(A.6) Quantum Sklyanin algebra. The fundamental commutation relation (A.1) is equivalent to the Sklyanin quadratic algebra which, in the general elliptic case, involves four generatorsŜ a , with a ∈ {1, 2, 3, 4}, and three 'structure constants' J ab ≡ −(J a − J b )/J c (with indices a, b, c mutually distinct) obeying the constraint J 12 + J 23 + J 31 + J 12 J 23 J 31 = 0. The Boltzmann weights specify an algebraic curve . The elliptic Sklyanin algebra admits two independent quadratic Casimir operatorsĈ (A.8) We are interested specifically in the trigonometric limit J 12 → 0, where the Sklyanin algebra reduces to the quantum algebra U q (su 2 ), representing a q-deformed quantum spin. In this limit, one findsĈ 0 −Ĉ 1 = 1, which leaves us with a single non-trivial Casimir operator of the formĈ the defining commutation relations of the trigonometric Sklyanin algebra read One can readily verify that these are indeed equivalent to standard U q (su 2 ) algebraic relations upon introducing the q-deformation parameter q = (1+κ)/(1−κ) and rescaling the generators, (A.13)

A.1 Semiclassical limit
We now perform the semiclassical limit. This is achieved, at the algebraic level, by expanding the fundamental commutation relation about the 'classical point' q = 1. Following Refs. [70,114], we parameterize q = e η , put η = and subsequently expand the quantum R-matrix to the leading order in , where 10 10 There is freedom to choose the overall scale of the r-matrix. Here we adopt the common convention [61], which amounts to fixing a particular normalization of the Sklyanin Poisson algebra, cf. Eq. (A.23). Notice that our normalization differs from that of Ref. [70].
is known as the classical r-matrix [70,75,114]. The classical amplitudes (Boltzmann weights) w a (λ) are similarly found by expanding the W -functions W a (λ) = i w a (λ) + O( 2 ), and structure constatnts J a = 1 − J a + O( 2 ), yielding The classical structure constants now obey J 12 + J 23 + J 31 = 0, and prescribe an algebraic where we have simultaneously introduced classical 'Sklyanin variables' 12 with classical structure constants J ab ≡ J a − J b . Poisson algebra (A.23) admits two independent Casimir functions By prescribing their values c 0 and c 1 , respectively, variables S a takes values on a twodimensional non-degenerate Poisson submanifold. By fixing them to the manifold is topologically equivalent to a two-sphere; {S a } 3 a=1 lie on a two-sphere, whereas S 0 > 0 is uniquely determined by S 3 through C 1 .
Introducing linear combinations K ± ≡ S 0 ± S 3 and putting K ≡ K + , we have K + K − = C 1 , yielding [115] A q : {K, S ± } = ∓i S ± K, The above Poisson relation also follow from quantum algebraic relations (A.11) upon rescaling κ as κ = , applying the canonical corespondence (A. 19), and finally sending → 0. Poisson algebra A q , with a pair of Casimir functions C 0 , C 1 , can be understood as the classical counterpart of U q (su 2 ).

B Factorization in Weyl variables
We provide explicit prescriptions for the elementary exchanges of Weyl matrix U and V that appear in the decomposition of the Lax matrix (cf. Section 2), written in the form of symplectic transformations on Weyl pairs (x , y ): and in 'Trotter parametrization' (2.41) 13 , (B.9) , (B.12) and in 'Trotter parametrization' (2.41) The combined U − V exchange (in Trotter parametrization (2.41)) gives the two-body map Φ

C Isospectral flows and local conservation laws
A key implication of the discrete zero-curvature condition (2.3) is the existence of infinitely many conserved quantities which are in mutual involution. As we shortly demonstrate, these are a corollary of isospectrality of the transfer function. We accordingly introduce a staggared monodromy matrix T τ (λ) ≡ T τ ({S }; λ), defined as a right-to-left ordered product of Lax matrices along a zig-zag path (i.e. initial sawtooth) of length L T τ (λ) = L S L ; λ + L S L−1 ; λ − · · · L S 2 ; λ + L S 1 ; λ − . (C.1) where λ ± ≡ λ±τ /2. Here and subsequently we leave the dependence on anisotropy implicit. By tracing out the auxiliary C 2 space, we obtain the transfer function t τ (λ) = Tr T τ (λ). (C.2) Isospectral evolution. By virtue of the zero-curvature condition (2.3) and cyclic property of the trace, it follows that which in combination imply invariance under the time propagation (2.6) The quadratic bracket (A.20) can be immediately lifted onto the level of monodromies By subsequently taking the trace over C 2 , one readily finds that transfer functions Poisson commute, t τ (λ), t τ (λ ) = 0, (C. 6) for any pair of complex spectral parameters λ, λ ∈ C and fixed τ . Preservation of the characteristic polynomial of T τ (λ) is referred to as isospectrality.
Local conservation laws. Because the transfer function t τ (λ) depends analytically on λ, it can be used as a generating function of integrals of motion. In the thermodynamic limit, L → ∞, this provides infinitely many functionally independent conserved quantities. Expanding t τ (λ) in a power series in λ nonetheless produces complex-valued and non-local conserved quantities. In order to construct real conservation laws, it is useful to utilize the following involutions of the Lax matrix L(λ) = σ 2 L(λ) σ 2 , L T (λ) = σ 2 L(−λ) σ 2 , (C.7) which likewise apply to the monodromy matrix T τ (λ). It follows that the moduli of the traces are also in mutual involution, and their expansion in λ thus yields manifestly real conserved quantities. One can in fact extract two seperate infinite towers of local conserved quantities, denoting them by Q ± k;τ , using the standard approach (see, for example, Ref. [61]) by expanding the logarithmic derivatives of the modulus of t τ (λ) around either of the two distinguished points λ ± 0 = i ∓ τ 2 at which the Lax matrix degenerates (becoming a projector of rank one), Det L λ ± 0 = 0, L λ ± 0 = |α ± β ± | , (C.8) yielding two independent sequences of local charges (C.9) Taking the continuous time limit τ → 0, both branches merge into the local conseration laws of the anisotropic lattice Landau-Lifshitz model.
Parity symmetry. The transfer function t τ (λ) is Z 2 -invariant under the spin-reversal transformation (π-rotation of the canonical spin S around the x-axis) implying F 1 (K) = K −1 . At the level of the Lax matrix, F 1 represents conjugation by Pauli matrix σ 1 , which implies invariance F 1 (t τ (λ)) = t τ (λ) and hence charges Q ± k;τ are all invariant under F 1 .

D Reductions
The anisotropic model of interacting spin in discrete space-time, introduced in Section 2, admits several important reductions. There are three distinct parameters that can be taken to zero: anisotropy , discrete time-step τ , and lattice spacing ∆ (which we have kept fixed ∆ = 1 throughout the text). The full set of models which descend from Φ τ, can be therefore arranged on the vertices of a cube, as depicted in Figure 10. We briefly describe them below. 14 generates a two-body symplectic propagator [56] Φ τ (S 1 , S 2 ) = 1 s 2 + τ 2 s 2 S 1 + τ 2 S 2 + τ S 1 × S 2 , s 2 S 2 + τ 2 S 1 + τ S 2 × S 1 ,