# Anisotropic Landau-Lifshitz model in discrete space-time

Žiga Krajnik, Enej Ilievski, Tomaz Prosen, Vincent Pasquier

► To cite this version:

Žiga Krajnik, Enej Ilievski, Tomaz Prosen, Vincent Pasquier. Anisotropic Landau-Lifshitz model in discrete space-time. SciPost Physics, 2021, 11 (3), pp.051. 10.21468/SciPostPhys.11.3.051 . cea-04525448

HAL Id: cea-04525448

https://cea.hal.science/cea-04525448

Submitted on 28 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Anisotropic Landau-Lifshitz model in discrete space-time

Žiga Krajnik¹*, Enej Ilievski¹, Tomaž Prosen¹ and Vincent Pasquier²

1 Faculty of Mathematics and Physics, University of Ljubljana,
Jadranska 19, 1000 Ljubljana, Slovenia

2 Institut de Physique Théorique, Université Paris Saclay, CEA,
CNRS UMR 3681, 91191 Gif-sur-Yvette, France

ziga.krajnik@fmf.uni-lj.si

## Abstract

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.

<figure><img src="image_2.png" /></figure>

Copyright Z. Krajnik *et al*.

This work is licensed under the Creative Commons

Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 05-05-2021

Accepted 25-08-2021

Published 09-09-2021

doi:10.21468/SciPostPhys.11.3.051

<figure><img src="image_3.png" /></figure>

Check for
updates

## Contents

1 **Introduction** 1

2 **Anisotropic symplectic spin model in discrete space-time** 2

2.1 Sklyanin Lax matrix 5

2.2 Factorization 7

2.3 Two-body propagator 9

3 **Magnetization transport** 11

3.1 Transport coefficients 12

3.2 Numerical simulations: easy-axis regime 13

3.3 Numerical simulations: easy-plane regime 16

4 **Conclusion** 18

**Appendices** 19

A **Semiclassical limit of quantum algebra** 19

B **Factorization in Weyl variables** 22
<table><tr><td>C</td><td>Isospectral flows and local conservation laws</td><td>23</td></tr><tr><td>D</td><td>Reductions</td><td>25</td></tr><tr><td colspan="2">References</td><td>26</td></tr></table>

# 1 Introduction

The advent of powerful computational tools [1, 2] has tremendously advanced our understanding of many-body physics in and out of equilibrium [3–6]. Dealing with thermodynamic systems of strongly interacting degrees of freedom nonetheless still presents a daunting 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–12].

Concurrently, we have seen an increased theoretical interest in various facets of nonequilibrium physics. Integrable systems have been in the spotlight lately [13–16, 16–22], largely due to their inherent non-ergodic features [23–26] and anomalous transport properties [5, 27–37] as recently covered in a compilation of review articles [38–43]. The study of nonequilibrium properties in classical integrable dynamical systems of interacting particles or fields has received comparatively less attention [44–51].

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 [52, 53] and classical cellular automata [42, 54] have recently attracted much interest. An integrable space-time discretization of the isotropic classical Heisenberg ferromagnet has been recently obtained in Ref. [55], and subsequently generalized to a large class of ‘matrix models’ [53] 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. [55], representing a ‘brick-wall type’ circuit composed from elementary two-body symplectic maps. Our model can be thus regarded as an ‘integrable Trotterization’ [53, 55–57] of the anisotropic lattice Landau-Lifshitz field theory [58–60] (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 constructed yet; it has only appeared previously in a paper by Hirota [61] in a bilinear form in terms of the discrete tau function. The outcome of our construction is an explicit finite-
step integrable integration scheme that can be efficiently employed to computationally study dynamical properties of the model. As an application, we compute the Drude weight and diffusion constant as functions of spin chemical potential (i.e. magnetization) and anisotropy.

**Outline.** This paper is structured as follows. We begin in Section 2 by describing the formal setting. First, we introduce a discrete zero-curvature property of an auxiliary linear transport problem on a (tilted, light-cone) space-time lattice and interpret it in terms of a local dynamical map acting on a pair of classical spins. We proceed by introducing a dynamical system on a discrete space-time lattice in the form of a 'symplectic circuit'. Next, in Section 2.1 we define the local phase space equipped with a Poisson structure, and introduce the Lax matrix of the lattice Landau-Lifshitz model. Moving on, in Section 2.2 we outline how to explicitly solve the zero-curvature relation by exploiting the underlying algebraic relations, yielding the local time-propagator of the model. Section 3 is devoted to a physics application, where we carry out a detailed numerical study of magnetization transport in grand-canonical Gibbs states for the easy-axis and easy-plane regimes of the model. In Section 4 we conclude with a brief summary of the main results.

## 2 Anisotropic symplectic spin model in discrete space-time

Before delving into the specifics of the model, we first introduce the general setting. We shall mostly follow the lines and use the notation from previous works on related discrete-time models with isotropic interactions [53, 55].

**Discrete space-time.** The physical space-time lattice is a two-dimensional square lattice with nodes $(\ell, t) \in \mathbb{Z}^2$, with $\ell$ and $t$ referring to the spatial and temporal indices, respectively. Throughout the paper we adopt the convention that time flows upwards while the spatial axis is oriented towards the right. Moreover, we adopt periodic boundary conditions in space by identifying $\ell \equiv \ell+L$ and additionally assume, for definiteness, the system length $L$ to be even.

Each site of the space-time lattice is attached a local physical degree of freedom $S_{\ell}^{\ell}$, which we consider to be a classical spin $S = (S^1, S^2, S^3)$ of unit length, $S \cdot S = S^+S^- + (S^3)^2 = 1$ (with $S^\pm \equiv S^1 \pm iS^2$) obeying the canonical ultralocal Lie-Poisson brackets

$$ \{S_{\ell}^{a}, S_{\ell'}^{b}\} = \delta_{\ell,\ell'} \sum_{c} \epsilon_{abc} S_{\ell}^{c}, \qquad (1) $$

where we have used the Levi-Civita symbol $\epsilon_{abc}$.

We furthermore introduce a square light-cone lattice by tilting the space-time lattice by 45° degrees, and assign to its vertices (nodes) $(n, m) \in \mathbb{Z}^2$ auxiliary variables $\phi_{n,m}$. Physical spins accordingly sit at the middle of the edges on the light-cone lattice.

**Discrete zero-curvature condition.** Integrability of a dynamical systems in discrete time is a direct manifestation of a discrete zero-curvature property of an associated linear transport problem [62, 63] for auxiliary variables $\phi$, in analogy to completely integrable Hamiltonian systems (see e.g. Refs. [60, 64, 65]). To realize parallel transport on the auxiliary light-cone lattice we introduce the lattice shift operators [66] along the light-cone directions (i.e. characteristics $\ell \pm t = \text{const}$)

$$ \phi_{n+1,m} = L_{n,m}^{(+)}(\lambda)\phi_{n,m}, \qquad \phi_{n,m+1} = L_{n,m}^{(-)}(\lambda)\phi_{n,m}. \qquad (2) $$

Here the local 'matrix propagators' $L^{(\pm)}$ represent certain matrix-valued functionals of physical variables which, in addition, depend analytically on a free complex parameter $\lambda$. In order to
<figure><img src="image_4.png" /><figcaption>Figure 1: Elementary plaquette of the discrete light-cone lattice: classical spins $S \in S^2$ (circles at the middle of the edges) are situated on the vertices of the discrete spacetime lattice. Primed variables $S'_{1,2}$ pertain to spins $S_{1,2}$ that are time-shifted forward by one unit of time by applying the propagator $\Phi$ (red square).</figcaption></figure>

ensure consistency of such parallel transport, there should be no ambiguity in the order of propagation when passing from $\phi_{n,m}$ to $\phi_{n+1,m+1}$. In this case, the shift matrices $L^{(\pm)}$ constitute the Lax pair and $\lambda$ 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^{\pm} \equiv L(\lambda^{\pm}; S)$, i.e. using a pair of spectral parameters $\lambda^{\pm}$ depending only on the light-cone direction (but not on position $\ell$ or time $t$).

Commutativity of the light-cone shifts is neatly encapsulated by the following discrete zero-curvature property [67, 68] around an elementary square plaquette of the light-cone lattice,<fnref n="1" />

$$L(\lambda^{+}; S_2)L(\lambda^{-}; S_1) = L(\lambda^{-}; S'_2)L(\lambda^{+}; S'_1), \qquad (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(\lambda^{\pm})$ 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 fulfilled 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 $\Phi : (S_1, S_2) \to (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 $\Phi$ defined implicitly via the zero-curvature property. In the following, we first outline how to use $\Phi$ as the elementary building block of a many-body dynamical map in the form of a 'brick-wall circuit'.

<footnote n="1">¹A magnetic field could be incorporated in the construction by adding constant 'twist' matrices to the zero-curvature condition. See [53] for a similar construction.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 2: Fabric of discrete space-time: the square space-time lattice comprising classical spin degrees of freedom $S_{\ell}^{t}$ (circles) superimposed on a square light-cone lattice represented by a tilted checkerboard. In the middle of each plaquette, there is a two-body symplectic map $\Phi$ (red square) which maps the pairs of physical spin variables forward in time.</figcaption></figure>

To this end, let $M_1 \cong 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 \equiv M_1^{\times L}$, and accordingly we introduce a many-body dynamical map $\Phi^{\text{full}}: M_L \to M_L$.

By virtue of the light-cone geometry, the full time-dynamics consists of an altering sequence of ‘odd’ and ‘even’ local propagators,

$$ (S_{2\ell-1}^{2t+2}, S_{2\ell}^{2t+2}) = \Phi(S_{2\ell-1}^{2t+1}, S_{2\ell}^{2t+1}), \quad (S_{2\ell}^{2t+1}, S_{2\ell+1}^{2t+1}) = \Phi(S_{2\ell}^{2t}, S_{2\ell+1}^{2t}), \quad (4) $$

respectively. This results in a staggered circuit as depicted in Figure 2. By making use of the embedding prescription,

$$ \Phi^{(j)} = I^{\times(j-1)} \times \Phi \times I^{\times(L-j-1)}, \quad (5) $$

for $j = 1, \dots, L-1$, where $I: M_1 \to M_1$ stands for a local (single-site) unit map $I(S) \equiv S$, ($\Phi^{(L)}$ correspondingly affects the $L$th and the first spin), the full propagator $\Phi^{\text{full}}$ for two units of time $t \mapsto t+2$ decomposes as

$$ \Phi^{\text{full}} = \Phi^{\text{even}} \circ \Phi^{\text{odd}}, \quad (6) $$

with odd and even propagators further factorizing into commuting two-body maps,

$$ \Phi^{\text{odd}} = \prod_{\ell=1}^{L/2} \Phi^{(2\ell)}, \quad \Phi^{\text{even}} = \prod_{\ell=1}^{L/2} \Phi^{(2\ell-1)}. \quad (7) $$

So far our construction has been entirely formal. In what follows, we derive the map $\Phi^{\text{full}}$ that corresponds to an integrable space-time discretization of the anisotropic lattice Landau-Lifshitz model, as introduced by Sklyanin [59, 60, 69]. In this view, $\Phi^{\text{full}}$ represents its ‘integrable Trotterization’.
## 2.1 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

$$ \mathcal{A}_q: \quad \{K, S^{\pm}\} = \mp i\rho S^{\pm}K, \quad \{S^+, S^-\} = -\frac{i\rho}{2}(K^2 - K^{-2}). \qquad (8) $$

This Poisson algebra admits a quadric (Casimir) function

$$ C_0 \equiv S^+S^- + \left( \frac{K-K^{-1}}{2} \right)^2, \qquad (9) $$

satisfying $\{C_0, S^a\} = 0$. In terms of ‘Sklyanin spins’ $^2$, with components $\{S^a\}_{a=0}^3$, with

$$ S^0 \equiv \frac{1}{2}(K + K^{-1}), \quad S^1 \equiv \frac{1}{2}(S^+ + S^-), \quad S^2 \equiv \frac{1}{2i}(S^+ - S^-), \quad S^3 \equiv \frac{1}{2}(K - K^{-1}), \qquad (10) $$

the Hamiltonian of the lattice Landau–Lifshitz model takes a compact form

$$ H_{\text{LLL}}^{\rho} \approx \sum_{\ell=1}^{L} \log \left[ \sum_{a=0}^{3} g_a S_{\ell}^{a} S_{\ell+1}^{a} \right], \qquad (11) $$

with couplings $g_0 = \sinh^2(\rho)$, $g_1 = g_2 = 1$ and $g_3 = \cosh^2(\rho)$.

We note that $\mathcal{A}_q$ is just the trigonometric limit of ‘elliptic spins’ that constitute the so-called quadratic Poisson–Sklyanin algebra, cf. Appendix A for further information. The algebra $\mathcal{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 $\rho \in \mathbb{R}_+$,

(ii) the easy-plane regime, with an interaction parameter $\gamma$ having a compact support $^3 \gamma \in [-\pi/2, \pi/2]$. The associated Hamiltonian is obtained from Eq. (11) by analytically continuing $\rho$ from the real axis onto the imaginary axis, $\rho \to i\gamma$.

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(\rho). \qquad (12) $$

In this way, we select a two-dimensional non-degenerate Poisson submanifold $M_{ea}$ inside $\mathbb{R}^3$ which is diffeomorphic to a 2-sphere. In the easy-plane regime, obtained under analytic continuation $\rho \to i\gamma$, the submanifold $M_{ep}$ is also topologically equivalent to a 2-sphere provided $\gamma \in [-\pi/2, \pi/2]$. The upshot here is that there exist a smooth bijective mapping to the symplectic sphere $S^2$, implying that variables $K$ and $S^\pm$ can be parameterized in terms of classical canonical spins $S$. A local mapping from $S^2$ to $M_{ea}$ is given by [60]

$$ \mathcal{K} = e^{\rho} S^3, \quad S^{\pm} = F_{\rho}(S^3)S^{\pm}, \qquad (13) $$

<footnote>$^2$Classical canonical spins can be recovered in the $\rho \to 0$ limit by simultaneously rescaling $S^a \to S^a/\rho$ for $a \in \{1, 2, 3\}$.</footnote>

<footnote>$^3$This restriction has been overlooked in the previous works [70, 71] on the lattice Landau–Lifshitz model in continuous time, which specialize to the value $\gamma = 1$.</footnote>
where

$$F_{\varrho}(s) \equiv \sqrt{\frac{\sinh^2(\varrho) - \sinh^2(\varrho s)}{1 - s^2}}, \qquad (14)$$

the form of which follows from the Casimir function, see Eqs. (9), (12). Two imporant remarks are in order at this stage: (I) in the *easy-plane regime* (with $\varrho \to i\gamma$ and $F_{i\gamma}(s)$), Eq. (13) provides a *bijective* map only inside a compact interval of anisotropies $|\gamma| \le \pi/2$; (II) there exist other Poisson submanifolds (classified in [60, 72]) 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 (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

$$L(z) = \frac{2}{z - z^{-1}} \begin{bmatrix} \frac{1}{2}(z K - z^{-1} K^{-1}) & S^{-} \\ S^{+} & -\frac{1}{2}(z^{-1} K - z K^{-1}) \end{bmatrix}, \qquad (15)$$

where we have for convenience introduced a multiplicative spectral parameter $z \in \mathbb{C}$,

$$z(\lambda) = e^{i\varrho\lambda}. \qquad (16)$$

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. (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 $\Phi$.<fnref n="4" /> This difficulty can be elegantly overcome by exploiting certain factorization properties of the Lax matrix into more fundamental constituents.

## 2.2 Factorization

In this section, we outline how the Lax matrix can be further decomposed into more elementary constituents. Such factorization properties are indeed deeply 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 $\mathcal{R}$-matrix admits a triangular (Borel) decomposition [73–76]. The resulting ‘partonic’ Lax operators can be used to construct the Baxter Q-operators [77]. There exist various equivalent constructions in the literature, see e.g. Refs. [78–80] and Refs. [81–83]. 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\varrho x y, \qquad (17)$$

<footnote n="4">⁴This is to some extent true even in the simplest isotropic case, where nonetheless it is possible to make use of vector identities to solve the zero-curvature condition explicitly.</footnote>
via the prescription

$$ \mathcal{K} = \nu y, \quad \mathcal{S}^+ = \frac{x}{2}(\nu^2 y - \nu^{-2}y^{-1}), \quad \mathcal{S}^- = -\frac{1}{2x}(y - y^{-1}), \qquad (18) $$

where we have simultaneously defined, for compactness of notation, $\nu \equiv e^{\mathcal{L}}$. We accordingly make the identification

$$ y = \nu^{-1} \mathcal{K} = e^{\mathcal{L}(S^3-1)}. \qquad (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 $\mathcal{S}^+ = \bar{\mathcal{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

$$ |x|^2 = \begin{cases} -[y]_{\nu^2}, & y \in \mathbb{R}_+ \\ +[y]_{\nu^2}, & y \in U(1), \end{cases} \quad [y]_{\nu^2} \equiv \frac{y - y^{-1}}{\nu^2 y - \nu^{-2}y^{-1}}. \qquad (20) $$

**Factorization.** By introducing<fnref n="5" />

$$ L_{x,y} \equiv \frac{1}{2}(z - z^{-1})L(z), \qquad (21) $$

along with a pair of spectral parameters,

$$ u = \nu z, \quad v = \nu z^{-1}, \qquad (22) $$

satisfying $uv = \nu^2$, $u/v = z^2$, the Lax operator from Eq. (15) can be factorized as

$$ L_{x,y}(V,U) \equiv \frac{1}{2} X^{-1} V Y U^T X = \frac{1}{2} \begin{bmatrix} u y - u^{-1} y^{-1} & -(y - y^{-1})x^{-1} \\ (u v y - u^{-1} v^{-1} y^{-1})x & -v y + v^{-1} y^{-1} \end{bmatrix}, \qquad (23) $$

in terms of diagonal ‘coordinate’ matrices $X$ and $Y$

$$ X = \begin{bmatrix} x & 0 \\ 0 & 1 \end{bmatrix}, \quad Y = \begin{bmatrix} y & 0 \\ 0 & y^{-1} \end{bmatrix}, \qquad (24) $$

and the ‘spectral’ matrices $U$ and $V$

$$ U = \begin{bmatrix} u & -u^{-1} \\ -1 & 1 \end{bmatrix}, \quad V = \begin{bmatrix} 1 & 1 \\ v & v^{-1} \end{bmatrix}. \qquad (25) $$

**Elementary exchange relations.** We next describe various elementary exchange procedures involving $L_{x,y}(U, V)$. Firstly, interchanging the spectral matrices

$$ L_{x,y}(V,U) = L_{x',y'}(U,V), \qquad (26) $$

specifies a map $(x, y) \mapsto (x', y')$. The canonical Weyl relations are preserved<fnref n="6" /> by additionally assuming $y' = y$, in which case

$$ X^{-1} V Y U^T X = X'^{-1} U Y V^T X'. \qquad (27) $$

<footnote n="5">⁵An overall multiplicative scalar factor bears no consequence on the dynamics as it cancels out from the discrete zero-curvature condition (3).</footnote>

<footnote n="6">⁶This can be inferred directly using that $x' = x f(y)$ and $y' = y$ preserve the Poisson bracket (17).</footnote>
Secondly, there exist quadratic exchange relations of the type

$$L_{x_1, y_1}(U_1, V_1)L_{x_2, y_2}(V_2, U_2) = L_{x'_1, y'_1}(U_1, V_2)L_{x'_2, y'_2}(V_1, U_2), \quad (28)$$

$$L_{x_1, y_1}(V_1, U_1)L_{x_2, y_2}(U_2, V_2) = L_{x'_1, y'_1}(V_1, U_2)L_{x'_2, y'_2}(U_1, V_2), \quad (29)$$

which are canonical provided that $x_1 = x'_1$, $x'_2 = x_2$, implying

$$Y_1 V_1^T X_1 X_2^{-1} V_2 Y_2 = Y_1' V_2^T X_1 X_2^{-1} V_1 Y_2', \quad (30)$$

$$Y_1 U_1^T X_1 X_2^{-1} U_2 Y_2 = Y_1' U_2^T X_1 X_2^{-1} U_1 Y_2'. \quad (31)$$

From Eq. (27) we can deduce the following transformation of $(x, y) \mapsto (x', y' = y)$:

$$x' = -x \frac{uvy - (uvy)^{-1}}{y - y^{-1}}. \quad (32)$$

Equations (30) and (31) are equivalent under the interchange $u \leftrightarrow v$, so it suffices to solve only one of them. By direct calculation, exchanging $V_1$ and $V_2$ according to Eq. (30) yields

$$(y_1')^2 = y_1^2 \frac{v_1}{v_2} \frac{v_1 x_2 + v_2 x_1}{v_1 x_1 + v_2 x_2}, \quad (y_2')^2 = y_2^2 \frac{v_2}{v_1} \frac{v_1 x_1 + v_2 x_2}{v_1 x_2 + v_2 x_1}, \quad x_\ell' = x_\ell, \quad (33)$$

$\ell \in \{1, 2\}$, while analogously the exchange of $U_1$ and $U_2$ according to Eq. (31) yields the same expression with $v_\ell$ replaced by $u_\ell$:

$$(y_1')^2 = y_1^2 \frac{u_1}{u_2} \frac{u_1 x_1 + u_2 x_2}{u_2 x_1 + u_1 x_2}, \quad (y_2')^2 = y_2^2 \frac{u_2}{u_1} \frac{u_2 x_1 + u_1 x_2}{u_1 x_1 + u_2 x_2}, \quad x_\ell' = x_\ell. \quad (34)$$

By the same reasoning as previously (only swapping the roles of $x_\ell$ and $y_\ell$), transformations (33) and (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

$$L_{x_1, y_1}(V_1, U_1)L_{x_2, y_2}(V_2, U_2) = L_{x'_1, y'_1}(V_1, U_2)L_{x'_2, y'_2}(V_2, U_1). \quad (35)$$

This can be achieved by a composition of three elementary exchanges

$$(V_1, U_1), (V_2, U_2) \rightarrow (V_1, U_1), (V_2, U_2) \rightarrow (V_1, U_2), (V_2, U_1). \quad (36)$$

In an analogous fashion we can exchange $V_1$ with $V_2$

$$L_{x_1, y_1}(V_1, U_1)L_{x_2, y_2}(V_2, U_2) = L_{x'_1, y'_1}(V_2, U_1)L_{x'_2, y'_2}(V_1, U_2), \quad (37)$$

via the following sequence of elementary exchanges

$$(V_1, U_1), (V_2, U_2) \rightarrow (U_1, V_1), (V_2, U_2) \rightarrow (U_1, V_2), (V_1, U_2) \rightarrow (V_2, U_1), (V_1, U_2). \quad (38)$$

**Solving the zero-curvature condition.** To obtain the two-body propagator $\Phi$ the discrete zero-curvature condition, Eq. (3),

$$L_{x_1, y_1}(V_1, U_1)L_{x_2, y_2}(V_2, U_2) = L_{x'_1, y'_1}(V_2, U_2)L_{x'_2, y'_2}(V_1, U_1), \quad (39)$$

can now be solved by composing a 'U-exchange' (35) with a 'V-exchange' (37), in whichever order

$$(V_1, U_1), (V_2, U_2) \rightarrow (V_1, U_2), (V_2, U_1) \rightarrow (V_2, U_2), (V_1, U_1). \quad (40)$$

Their composition gives an explicit propagator $\Phi$ that solves the zero-curvature condition (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.
## 2.3 Two-body propagator

As in the isotropic case [55], the discrete 'Trotter' time-step $\tau \in \mathbb{R}$ enters the construction through the difference of the additive alternating (i.e. staggered) spectral parameters $\lambda^{\pm} \equiv \lambda \pm \tau/2$, carried by the light-cone Lax matrices, cf. Eq. (3). In terms of spectral parameter $u$ and $v$, this amounts to setting

$$u^{\pm} = v e^{i\varrho} \lambda^{\pm}, \quad v^{\pm} = v e^{-i\varrho} \lambda^{\pm}, \qquad (41)$$

recalling that $v = e^{i\varrho}$. With the above prescription, the exchange relations (given by equations (35) and (37)) and the propagator $\Phi$ simplify considerably. Explicit expressions can be found in Appendix B).

The elementary propagator $\Phi_{\tau,\varrho}$ can alternatively be written explicitly using Sklyanin variables $K, S^{\pm}$. With aid of Eqs. (18) and some exercise, we find (where $w = e^{i\varrho} \tau$)

$$\begin{gathered} \left(\frac{K'_1}{K_1}\right)^2 = \frac{(v^2 + v^{-2})(K_1 K_2 + K_1^{-1} K_2^{-1}) - 2K_1 K_2^{-1} - (w^2 + w^{-2})K_1^{-1} K_2 + 4w S_1^+ S_2^- + 4w^{-1} S_1^- S_2^+}{(v^2 + v^{-2})(K_1 K_2 + K_1^{-1} K_2^{-1}) - 2K_1^{-1} K_2 - (w^2 + w^{-2})K_1 K_2^{-1} + 4w S_1^- S_2^+ + 4w^{-1} S_1^+ S_2^-}, \\ \left(\frac{K'_2}{K_2}\right)^2 = \frac{(v^2 + v^{-2})(K_1 K_2 + K_1^{-1} K_2^{-1}) - 2K_1^{-1} K_2 - (w^2 + w^{-2})K_1 K_2^{-1} + 4w S_1^- S_2^+ + 4w^{-1} S_1^+ S_2^-}{(v^2 + v^{-2})(K_1 K_2 + K_1^{-1} K_2^{-1}) - 2K_1 K_2^{-1} - (w^2 + w^{-2})K_1^{-1} K_2 + 4w S_1^+ S_2^- + 4w^{-1} S_1^- S_2^+}, \\ \left[\frac{(S_1^{\pm})'}{(S_2^{\pm})'}\right] = \Omega(w^{\pm 1}) \left[\frac{S_1^{\pm}}{S_2^{\pm}}\right], \end{gathered} \qquad (42)$$

with

$$\Omega(w) = \frac{1}{K_1 K_2 w - (K_1 K_2 w)^{-1}} \left[ \frac{K'_1 K_2 - (K'_1 K_2)^{-1}}{(K'_2/K_2)w - (K_2/K'_2)w^{-1}} \frac{(K'_1/K_1)w - (K_1/K'_1)w^{-1}}{K'_1 K_2 - (K'_1 K_2)^{-1}} \right]. \quad (43)$$

We remind the reader that in the easy-axis regime $K \in \mathbb{R}_+$, whereas in the easy-plane regime $\gamma S^3 = \arg(K) \in [-\pi/2, \pi/2]$. Therefore, there is no ambiguity in Eqs. (42) when taking the square roots to determine $K'_\ell$.

**Properties.** There are a number of important properties worth highlighting:

* In the continuous-time limit $\tau \to 0$, the propagator $\Phi_{\tau,\varrho}$ reduces to the identity map.

* The mapping (42) manifestly preserves the product $K_1 K_2$, implying that the local propagator conserves the third component of the spin,

$$S_1^3 + S_2^3 = (S_1^3)' + (S_2^3)' \qquad (44)$$

This local $U(1)$ symmetry directly implies conservation of $S_{tot}^3 \equiv \sum_{\ell=1}^L S_\ell^3$ under the full propagator $\Phi^{\text{full}}$ which, in the Hamiltonian limit, $\tau \to 0$, becomes the Noether charge of the uniaxially anisotropic Landau-Lifshitz model.

* In the easy-axis regime, the local propagator (42) is a periodic function of interaction $\varrho$, with period $2\pi/\tau$.

* The two-body propagator $\Phi_{\tau,\varrho}$ is invariant under the transformation

$$\varrho \to -\varrho, \quad \tau \to -\tau, \quad \mathbf{S}_{1,2} \to -\mathbf{S}_{1,2}. \qquad (45)$$

It is therefore sufficient to consider only the interval $\varrho \in [0, \pi/\tau]$.

* In the easy-axis regime, at $\varrho = \pi/\tau$ (equiv. $w = -1$) the dynamics trivializes to the identity for all $\tau$,

$$\lim_{\varrho \to \pi/\tau} \Phi_{\tau, \varrho}(\mathbf{S}_1, \mathbf{S}_2) = (\mathbf{S}_1, \mathbf{S}_2). \qquad (46)$$
* In the easy-plane regime, the dynamics is not periodic in $\rho$, but is instead confined to the interval $\gamma \in [-\pi/2, \pi/2]$. By virtue of Eqs. (45), it is sufficient to consider only the range $\gamma \in [0, \pi/2]$.

* In the easy-plane regime, in the limit $\tau \to \infty$ with $\gamma$ fixed (i.e. $|w| \to \infty$), $\Phi_{\tau,i\gamma}$ reduces to a permutation

$$ \lim_{\tau \to \infty} \Phi_{\tau,i\gamma}(S_1, S_2) = (S_2, S_1). \qquad (47) $$

* The two-body propagator $\Phi$ is a classical Yang-Baxter map [66, 84-86], obeying the braided Yang-Baxter relation on the Cartesian product $M_1^{\times 3}$ (suppressing the anisotropy parameter)

$$ \Phi_{\lambda_2-\lambda_3}^{(1,2)} \circ \Phi_{\lambda_1-\lambda_3}^{(2,3)} \circ \Phi_{\lambda_1-\lambda_2}^{(1,2)} = \Phi_{\lambda_1-\lambda_2}^{(2,3)} \circ \Phi_{\lambda_1-\lambda_3}^{(1,2)} \circ \Phi_{\lambda_2-\lambda_3}^{(2,3)}, \qquad (48) $$

where $\Phi^{(i,j)}$ signifies that $\Phi$ operates non-identically on i-th and j-th copy of $M_1^{\times 3}$.

**Discrete Noether current.** The lattice Landau-Lifshitz model, cf. Eq. (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. (44), that the same property holds even in the discrete-time setting.

<figure><img src="image_6.png" /><figcaption>Figure 3: Discrete Noether current, shown at an even site $2\ell$: the time-increment $t \to t+2$ of charge density, $(q_{2\ell}^{t+2} - q_{2\ell}^t)/\tau$ equals the difference of two nearby current densities at consecutive times, $j_{2\ell+1}^{2t+1} - j_{2\ell}^{2t}$. Current densities $j_\ell^t$, depicted by pink arrows, can be compactly expressed as the difference of charge densities at consecutive times, $j_\ell^t = (q_\ell^{t+1} - q_\ell^t)/\tau$.</figcaption></figure>

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 $\Phi^{\text{full}}$, the discrete continuity equations at even and odd sites are of the form

$$ \frac{1}{\tau} (q_{2\ell}^{2t+2} - q_{2\ell}^{2t}) + j_{2\ell+1}^{2t+1} - j_{2\ell}^{2t} = 0, \qquad \frac{1}{\tau} (q_{2\ell+1}^{2t+2} - q_{2\ell+1}^{2t}) + j_{2\ell+2}^{2t} - j_{2\ell+1}^{2t+1} = 0, \qquad (49) $$

respectively. While the charge densities $q_\ell^t \equiv q(\ell, t)$ are ultralocal, namely they sit at site $\ell$ and
time $t$, the associated current density corresponds to their forward differences <fnref n="7" />

$$j_{\ell}^{t} \equiv j(\ell, t) = \frac{1}{\tau}(q(\ell, t+1) - q(\ell, t)), \qquad (50)$$

which can be immediately verified by plugging the expression into the discrete continuity equations (49) and using the conservation property (44). Notice that $q_{\ell}^{t+1}$ is a function of two adjacent charge densities at the previous time slice, namely $q_{\ell}^{t}$ and $q_{\ell \pm 1}^{t}$, cf. Figure 3, that enter the local propagator $\Phi$.

**Isotropic, continuous time and field-theory limits.** The two-parametric dynamical map $\Phi_{\tau, \varrho}^{\text{full}}$ admits various important limits: (i) the limit of vanishing anisotropy $\varrho \to 0$, (ii) the continuous-time limit $\tau \to 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’.

## 3 Magnetization transport

We demonstrate the utility of the constructed dynamical system by studying magnetization transport in thermal equilibrium. We focus on a 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.

### 3.1 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 $\rho_1(\mu)$ on $M_1$,

$$\rho_1(\mathbf{S}) = \frac{1}{4\pi \sinh \mu} e^{\mu S^3}, \qquad (51)$$

which can be interpreted as the grand-canonical measure with chemical potential $\mu$.<fnref n="8" /> We accordingly define a separable stationary measure $\rho_L \equiv \rho_L(\{\mathbf{S}_\ell\}_{\ell=1}^L)$ on the full phase space $M_L$,

$$\rho_L = \prod_{\ell=1}^{L} \rho_1(\mathbf{S}_\ell), \qquad (52)$$

with chemical potential $\mu$ enforcing a non-zero average of magnetization. Invariance of $\rho_L$ under the full dynamics $\Phi^{\text{full}}$ follows from the fact that the quadratic Sklyanin bracket is preserved by the action of $\Phi$ (see Ref. [53] for an analogous formal proof).

We now pass over to the thermodynamic limit by sending $L \to \infty$. The grand-canonical free energy per site is given by

$$f(\mu) = -\log Z_1(\mu), \qquad (53)$$

<footnote n="7">⁷We note that the current density is only defined up to a total difference. We fix this gauge freedom by demanding that the current vanishes when both input spins are equal.</footnote>

<footnote n="8">⁸Such a one-body measure maximizes the Shannon (or Gibbs) entropy of a state, that is $-\int d\Omega \rho_1(\mathbf{S}) \log \rho_1(\mathbf{S})$ where $d\Omega$ is a uniform measure on $S^2$ with a prescribed average $\langle S^3 \rangle$ set by a Lagrange multiplier $\mu$. In the continuous-time (i.e. Hamiltonian) limit, $\tau \to 0$, $\rho_1(\mathbf{S})$ is the grand-canonical Gibbs measure in the limit of infinite temperature.</footnote>
where

$$Z_1(\mu) = \int_{S^2} d\Omega \rho_1(\mathbf{S}) = \frac{4\pi \sinh \mu}{\mu}, \qquad (54)$$

represents the onsite grand-canonical partition function, with $d\Omega$ denoting the volume element on $S^2$. Writing $d\Omega^{\text{full}} = \prod_{\ell=1}^{L} d\Omega(S_\ell)$, the corresponding thermodynamic average of any local observable $\mathcal{O}$ is therefore computed as

$$\langle \mathcal{O} \rangle_{\mu} = \lim_{L \to \infty} [Z_1^L(\mu)]^{-1} \int_{M_L} d\Omega^{\text{full}} \rho_L \mathcal{O}. \qquad (55)$$

For example, the average of magnetization density reads

$$\langle S^3 \rangle_{\mu} = \coth \mu - 1/\mu. \qquad (56)$$

In the following we accordingly take $L \to \infty$ 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 $\sigma(\omega)$, which is in general decomposed as

$$\sigma(\omega) = \pi D \delta(\omega) + \sigma_{\text{reg}}(\omega). \qquad (57)$$

The $\delta$-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 $\sigma_{\text{reg}}(0)$, characterizes transport on sub-ballistic scales. Provided $\sigma_{\text{reg}}(0)$ is finite, it can be associated to the spin diffusion constant defined as

$$D = \chi^{-1} \lim_{\omega \to 0} \sigma_{\text{reg}}(\omega), \qquad (58)$$

where $\chi$ denotes static spin susceptibility, $\chi(\mu) = \partial^2 f(\mu)/\partial^2 \mu = 1 + \mu^{-2} - \coth^2 \mu$. In integrable interacting systems, the diffusion constant quantifies the diffusive spreading which quasiparticle excitations exhibit upon elastic collisions (see e.g. Refs. [40,87]). As customary, we define the dynamical structure factor

$$S(l, t) = \langle q(l, t) q(0, 0) \rangle_{\mu}^{c}, \qquad (59)$$

with a sum rule $\sum_{\ell \in \mathbb{Z}} S(l, t) = \chi$. The spin Drude weight governs the time-asymptotic growth of its second moment

$$\sum_{\ell \in \mathbb{Z}} \ell^2 S(l, t) \approx D t^2 + \mathcal{O}(t), \qquad (60)$$

where the subleading correction, which equals $2\chi D$, stores information about the spin diffusion constant $D$. We note that numerically extracting $D$ from the $\mathcal{O}(t)$ broadening of the second moment (60) 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 (49) and is given by

$$J(2t) = \sum_{\ell \in \mathbb{Z}} (j(2\ell, 2t) + j(2\ell + 1, 2t + 1)). \qquad (61)$$

The spin Drude weight corresponds to the asymptotic value of the current autocorrelation function

$$\mathcal{D} = \lim_{t \to \infty} \frac{1}{t} \sum_{t'=0}^{t} \langle J(2t')j(0, 0) \rangle_{\mu}^{c}, \qquad (62)$$
whereas the spin diffusion constant is the integrated current autocorrelator with the asymptotic value subtracted

$$D = \chi^{-1} \lim_{t \to \infty} \tau \sum_{t'=0}^{t} \left( \langle J(2t')j(0,0) \rangle_{\mu}^{c} - \mathcal{D} \right). \quad (63)$$

Lastly, we introduce the dynamical exponent $\alpha$ that characterizes the late-time decay of the dynamical spin structure factor

$$S(0, t) \sim t^{-\alpha}. \quad (64)$$

Ballistic spreading corresponds to exponent $\alpha = 1$. On the other hand, we can speak of normal diffusion when (i) $\alpha = 1/2$ and (ii) the dynamical structure factor converges asymptotically to a Gaussian stationary scaling profile,

$$S(l, t) \approx \frac{\chi}{\sqrt{2\pi Dt}} e^{-l^2/2Dt}. \quad (65)$$

<figure><img src="image_7.png" /><figcaption>Figure 4: (a) Spin Drude weight $\mathcal{D}$ and (b) diffusion constant $D$ in the discrete Landau-Lifshitz model, cf. Eq. 4, in the easy-axis regime, depending on the chemical potential $\mu$ and the anisotropy $\rho$. (a, inset) Drude weight (purple) and diffusion constant (olive green) as functions of $\mu$ at $\rho = 1$ (vertical grey line). (b, inset) Diffusion of $D$ at zero average magnetization, $\mu = 0$ (blue curve), in the isotropic limit $\rho \to 0$, compared to $\sim 1/\rho$ power-law scaling (black dashed line). Simulation parameters: $\tau = 1$, $L = 2^{12}$, averaged over $10^5$ initial configurations.</figcaption></figure>

## 3.2 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 $\rho$ and arbitrary finite magnetization density $\mu > 0$. In the unmagnetized ensemble ($\mu = 0$), corresponding to a uniform measure $\rho_1 = 1/\text{Vol}(S^2) = 1/4\pi$, the Drude weight vanishes. This point is distinguished by the discrete $\mathbb{Z}_2$ symmetry, corresponding to the global reversal of all spins, $S_\ell^\pm \mapsto S_\ell^\mp$, $S_\ell^z \mapsto -S_\ell^z$, which is a canonical transformation, i.e. preserving the Poisson bracket. We find, moreover, that both transport coefficients vanish when $\rho \to \pi/\tau$, consistent with the dynamics becoming trivial in that limit, cf. Eq. (46).
<figure><img src="image_8.png" /><figcaption>Figure 5: (a) Finite-time dynamical exponents $\alpha$ ($t \in [100, 300]$ dashed-dotted lines, $t \in [300, 1000]$ solid lines) depending on the anisotropy $\rho$ at four values of the chemical potentials $\mu$. Magnetized ensembles ($\mu > 0$) exhibit ballistic transport with exponent $\alpha = 1$. The sub-leading diffusive correction increases with $\rho$ and overwhelms the ballistic contribution at finite times. In the non-magnetized sector at $\mu = 0$, the dynamical exponent is approximately $\alpha \approx 2/3$ (black circle) when approaching $\rho \to 0$, consistent with KPZ-type superdiffusion at the isotropic point $\rho = 0$. At large values of $\rho$, the dynamical exponent approaches $\alpha \approx 1/2$ which is indicative of normal diffusion. (b) Rescaled diffusive stationary scaling profiles of the dynamical structure factors $t^{1/2}C(l, t)$ at $\mu = 0$ and different values of $\rho$. Simulation parameters: $\tau = 1, L = 2^{12}$, averaged over $10^5$ initial configurations.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Ballistic stationary scaling profiles of dynamical structure factors $S(l, t)$ (shown in insets) computed at $\mu = 2, 5$. Simulation parameters: $\tau = 1, L = 2^{12}$, averaged over $10^5$ initial configurations.</figcaption></figure>

Fixing anisotropy to $\rho = 1$, we find a non-monotonic dependence of $D$ on the chemical potential $\mu$, see the inset plot in Figure 4, which is consistent with the picture of dressed quasiparticles [40, 41]: in an unpolarized ensemble, the dressed magnetization vanishes and so does the spin Drude weight. For small $\mu$, the quasiparticles behave paramagnetically, with dressed magnetization behaving as $\sim \mu$. For strong polarizations, $\mu \to \infty$, 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 $\mu > 0$. We find that it decreases monotonically with increasing $\mu$. At any given anisotropy $\rho$, $D$ grows upon lowering $\mu$, and attains its maximal value at $\mu = 0$. There, upon approaching the isotropic point, $\rho \to 0$, the diffusion constant $D$ diverges in (approximately) algebraic fashion, $D \sim 1/\rho^{\chi_{ea}}$, as shown in Figure 4 (right inset). This conforms with a theoretical expectation $\chi_{ea} = 1$ based on the known behavior in the quantum Heisenberg chain [34].<fnref n="9" /> A reliable estimation of $\chi_{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 [88] which has recently been under intense scrutiny [30, 32, 53, 55, 71, 89–91]. There exists ample numerical evidence [32, 33, 53, 55, 71, 92, 93] that this phenomenon belongs to the universality class of the Kardar–Parisi–Zhang (KPZ) equation [94]. In Figure 5, we numerically corroborate the predicted anomalous dynamical exponent $\alpha = 2/3$.

<figure><img src="image_10.png" /><figcaption>Figure 7: (a) Spin Drude weight $D$ and (b) and diffusion constant $D$ in the discrete Landau–Lifshitz model, cf. Eq. (4), in the easy plane-regime, depending on chemical potential $\mu$ and anisotropy $\gamma$. (a, inset) Drude weight (purple) and diffusion constant (olive green) as functions of $\mu$ at $\gamma = 1$ (vertical grey line). (b, inset) Divergence of $D$ at zero average magnetization, $\mu = 0$ (blue curve), in the isotropic limit $\gamma \to 0$, compared to $\sim 1/\rho$ power-law scaling (black dashed line). Simulation parameters: $\tau = 1, L = 2^{12}$, averaged over $10^5$ initial configurations.</figcaption></figure>

The transient effects can be discerned by plotting the dynamical structure factors (the equal-time profiles of dynamical correlation functions of the spin density), as shown in Figure 5. While the dynamics is ballistic for all values of $\rho > 0$ and $\mu > 0$, the finite-time dynamical exponent estimated based on (64) is found to be ballistic, $\alpha = 1$, only below a threshold value of $\rho$ (which increases with increasing $\mu$). This is due to the fact that the value of $D$ is small in comparison with $D$. By contrast, at $\mu = 0$, the estimated exponent is approximately diffusive, $\alpha \approx 1/2$ for $\rho \gtrsim 1$, while when approaching $\rho \to 0$ it starts to drift towards the expected superdiffusive exponent $\alpha = 2/3$. At $\mu = 0$, the stationary cross-sections shown in Figure 5 (b) are single-peaked, consistent with the absence of ballistic transport. As $\rho \to 0$, the broadening of the central peak is compatible with the aforementioned divergence of the diffusion constant.

<footnote n="9">⁹In the gapped phase of the Heisenberg spin-1/2 chain with anisotropy $\Delta = \cosh(\eta)$, the spin diffusion constant diverges as $\sim 1/\eta$ as $\eta \to 0$. In the semiclassical limit, easy-axis anisotropy $\rho$ is proportional to $\eta$, see Appendix A.</footnote>
<figure><img src="image_11.png" /></figure>

Figure 8: (a) Finite-time dynamical exponents $\alpha$ (full blue line for $t \in [300, 1000]$, blue dashed-dotted line for $t \in [100, 300]$) as functions of anisotropy for different values of chemical potentials $\mu$. At all $\mu > 0$, we find ballistic transport ($\alpha = 1$). At $\mu = 0$ and $\gamma \to 0$, the dynamical exponent crosses over to superdiffusion with $\alpha = 2/3$ (black circle). (b) Rescaled ballistic stationary profiles of the dynamical structure factor $tC(l, t)$ at $\mu = 0$.

### 3.3 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 $\mu$ and $\gamma$, except at the isotropic point at zero magnetization ($\mu = \gamma = 0$). Upon lowering $\gamma$, the spin Drude weight $D$ now monotonously decreases, for any fixed $\mu$, to a limiting value at $\gamma = 0$ which is non-zero for $\mu \neq 0$.

<figure><img src="image_12.png" /></figure>

Figure 9: Rescaled ballistic stationary profiles of the dynamical structure factor $tC(l, t)$ for (a) $\mu = 2$ and (b) $\mu = 5$, with the corresponding space-time dependence $C(l, t)$ (inset plots). Simulation parameters: $\tau = 1$, $L = 2^{12}$, averaged over $10^5$ initial configurations.

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 nu-
merical studies of the lattice Landau–Lifshitz model [70, 71], 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 $\mu \to 0$ limit, mirroring the situation in the easy-plane phase of the quantum Heisenberg XXZ spin chain [41, 95],

(b) from the viewpoint of the Mazur–Suzuki hydrodynamic projection [96, 97] (for more information see e.g. [87, 98, 99]) a finite value of $\mathcal{D}$ necessitates the existence of local (or quasilocal, see Refs. [98, 100, 101]) conserved quantities with finite spin-charge susceptibility even at $\mu \to 0$, that is $\langle JQ \rangle_0^c > 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. [102–105] and [101] for a review). They are commonly referred to as the ‘Z-charges’. Presently however, the standard local conservation laws $Q_{k;\tau}^{\pm}$ extracted by expanding the trace of the monodromy matrix as an analytic series (see the derivation in Appendix C) all satisfy $\langle JQ_{k;\tau}^{\pm} \rangle_0^c = 0$. This indeed follows directly from a simple symmetry argument: while all local charges $Q_{k;\tau}$ 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 $\gamma$ provided $\mu > 0$. In the $\mu \to 0$ limit however, it once again diverges as $D \sim 1/\gamma^{x_{ep}}$ algebraically with an estimated exponent $x_{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 $\tau$ (with fixed $\gamma$) the two-particle map reduces to the permutation, see Eq. (47), meaning that the dynamics becomes non-interacting and hence the structure factor concentrates near the light-cone boundary peaked at $\ell/t = \pm 2$.

**4 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 per-
ceived 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
zero-curvature condition on an auxiliary light-cone lattice, which implicitly defines the ele-
mentary 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 pre-
serves integrability. As an important immediate physics application, we numerically studied
magnetization transport in grand-canonical stationary ensembles at different values of magne-
tization 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 [41] (characterized by the dynamical exponent $\alpha = 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. [70, 71]. It is well-known indeed that the same behavior occurs in 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 [101–103] 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 [106, 107] and diffusion constants [108–110] obtained within the formalism of generalized hydrodynamics [27, 28]. This computation however requires the knowledge of the quasiparticle content of the model and the associated thermodynamic state functions [87], which are yet to be derived, e.g. by means of the inverse scattering techniques [60, 64, 65, 111] adapted to the fully discrete setting. We leave this task for future research.

## Acknowledgements

The work has been supported by ERC Advanced grant 694544 – OMNES and the program P1-0402 of Slovenian Research Agency.

## 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 consider the so-called fundamental commutation (or RLL) relation of the difference form [69, 73, 112]

$$R_{12}(\lambda - \lambda')\hat{L}_{13}(\lambda)\hat{L}_{23}(\lambda') = \hat{L}_{23}(\lambda')\hat{L}_{13}(\lambda)R_{12}(\lambda - \lambda'), \quad (A.1)$$

representing an operatorial identity on the tensor product of three vector spaces of the form $V_{1/2} \otimes V_{1/2} \otimes 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 \otimes V_S$, where $V_S$ pertains to a unitary irreducible representation of ‘Sklyanin quantum spins’ $\hat{S}^a$ (with $a \in \{0, 1, 2, 3\}$) of dimension $2S + 1$. The Lax operator
assumes the form

$$L(u) = 1 \otimes \hat{S}^0 + \sum_{a=1}^{3} W_a(\lambda) \sigma^a \otimes \hat{S}^a. \qquad (A.2)$$

**Quantum R-matrix.** The R-matrix acts on two copies of $V_{1/2} \cong \mathbb{C}^2$ which are regarded as auxiliary fundamental spins (subsequently represented by Pauli matrices $\sigma^a$). The R-matrix swaps (intertwines) the spectral parameters of the Lax operators. The R-matrix has the 6-vertex form [112]

$$R(\lambda) = 1 + \sum_{a=1}^{3} W_a(\lambda) \sigma^a \otimes \sigma^a, \qquad (A.3)$$

with trigonometric amplitudes (Boltzmann weights)

$$W_1(\lambda) = W_2(\lambda) = \frac{\sinh(\eta)}{\sinh(\lambda + \eta)}, \qquad W_3(\lambda) = \frac{\tanh(\eta)}{\tanh(\lambda + \eta)}, \qquad (A.4)$$

where

$$q = e^{\eta} \in \mathbb{R}_{+} \qquad (A.5)$$

is the quantum deformation parameter. The 6-vertex quantum R-matrix, cf. Eq. (A.3), is a specialization of more general 8-vertex elliptic R-matrix introduced by Sklyanin [72], obtained by sending the elliptic modulus to zero.

When embedded in $V_{1/2} \otimes V_{1/2} \otimes V_{1/2}$, the R-matrix satisfies the quantum Yang-Baxter relation [77, 112]

$$R_{12}(\lambda - \lambda') R_{13}(\lambda) R_{23}(\lambda') = R_{23}(\lambda') R_{13}(\lambda) R_{12}(\lambda - \lambda'). \qquad (A.6)$$

**Quantum Sklyanin algebra.** The fundamental commutation relation (A.1) is equivalent to the Sklyanin quadratic algebra

$$[\hat{S}^0, \hat{S}^a] = i \mathcal{J}_{bc} [\hat{S}^b, \hat{S}^c]_+, \qquad [\hat{S}^a, \hat{S}^b] = i [\hat{S}^0, \hat{S}^c]_+, \qquad (A.7)$$

which, in the general elliptic case, involves four generators $\hat{S}^a$, with $a \in \{1, 2, 3, 4\}$, and three ‘structure constants’ $\mathcal{J}_{ab} = -(\mathcal{J}_a - \mathcal{J}_b)/\mathcal{J}_c$ (with indices $a, b, c$ mutually distinct) obeying the constraint $\mathcal{J}_{12} + \mathcal{J}_{23} + \mathcal{J}_{31} + \mathcal{J}_{12}\mathcal{J}_{23}\mathcal{J}_{31} = 0$. The Boltzmann weights specify an algebraic curve $W_a^2 - W_b^2 = \mathcal{J}_{ab}(W_c^2 - 1)$. The elliptic Sklyanin algebra admits two independent quadratic Casimir operators

$$\hat{C}_0 = \sum_{a=0}^{3} (\hat{S}^a)^2, \qquad \hat{C}_1 = \sum_{a=1}^{3} \mathcal{J}_a (\hat{S}^a)^2. \qquad (A.8)$$

We are interested specifically in the trigonometric limit $\mathcal{J}_{12} \to 0$, where the Sklyanin algebra reduces to the quantum algebra $\mathcal{U}_q(\mathcal{S}u_2)$, representing a q-deformed quantum spin. In this limit, one finds $\hat{C}_0 - \hat{C}_1 = 1$, which leaves us with a single non-trivial Casimir operator of the form

$$\hat{C}_q \equiv q^{-1} \hat{K}^2 + q \hat{K}^{-2} + (q - q^{-1}) \hat{S}^+ \hat{S}^- \qquad (A.9)$$

Moreover, writing $\kappa \equiv \sqrt{\mathcal{J}_{23}} = \sqrt{(\mathcal{J}_3 - \mathcal{J}_2)/\mathcal{J}_1}$ and introducing

$$\hat{S}^{\pm} \equiv \hat{S}^{1} \pm i \hat{S}^{2}, \qquad \hat{K} \equiv \hat{S}^{0} + \kappa \hat{S}^{3} = q^{\hat{S}^{3}}, \qquad (A.10)$$

the defining commutation relations of the trigonometric Sklyanin algebra read

$$[\hat{K}, \hat{S}^{\pm}] = \pm \kappa [\hat{K}, \hat{S}^{\pm}]_{+}, \qquad [\hat{S}^{+}, \hat{S}^{-}] = \frac{1}{\kappa}(\hat{K}^{2} - \hat{K}^{-2}). \qquad (A.11)$$
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+\kappa)/(1-\kappa)$ and rescaling the generators,

$$ \mathcal{S}^{\pm} \equiv \frac{1}{2} \sqrt{(1+\kappa)(1-\kappa)} \mathcal{S}^{\pm}, \qquad \mathcal{K} \equiv \mathcal{K}, \qquad (A.12) $$

yielding

$$ \mathcal{K} \mathcal{S}^{\pm} = q^{\pm 1} \mathcal{S}^{\pm} \mathcal{K}, \qquad [\mathcal{S}^{+}, \mathcal{S}^{-}] = \frac{\mathcal{K}^2 - \mathcal{K}^{-2}}{q - q^{-1}}. \qquad (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. [69, 113], we parameterize $q = e^{\eta}$, put $\eta = \hbar \rho$ and subsequently expand the quantum $R$-matrix to the leading order in $\hbar$,

$$ R(\lambda) = 1 - 2i\hbar r(\lambda) + \mathcal{O}(\hbar^2), \qquad (A.14) $$

where $^{10}$

$$ r(\lambda) = -\frac{1}{2} \sum_{a=1}^{3} w_a(\lambda) \sigma^a \otimes \sigma^a, \qquad (A.15) $$

is known as the classical $r$-matrix [69, 74, 113]. The classical amplitudes (Boltzmann weights) $w_a(\lambda)$ are similarly found by expanding the $W$-functions $W_a(\lambda) = i\hbar w_a(\lambda) + \mathcal{O}(\hbar^2)$, and structure constants $J_a = 1 - \hbar J_a + \mathcal{O}(\hbar^2)$, yielding

$$ w_{1,2}(\lambda) = \frac{\rho}{\sin(\lambda)}, \qquad w_3(\lambda) = \frac{\rho}{\tan(\lambda)}, \qquad (A.16) $$

with anisotropy parameter

$$ \rho \equiv \sqrt{J_3 - J_1} \in \mathbb{R}_+. \qquad (A.17) $$

The classical structure constants now obey $J_{12}+J_{23}+J_{31} = 0$, and prescribe an algebraic curve $w_a^2(\lambda) - w_b^2(\lambda) = J_b - J_a$.

The classical $r$-matrix satisfies the classical analogue of the Yang-Baxter relation [69, 74]

$$ [r_{12}(\lambda - \lambda'), r_{13}(\lambda)] + [r_{12}(\lambda - \lambda'), r_{23}(\lambda')] + [r_{13}(\lambda), r_{23}(\lambda')] = 0, \qquad (A.18) $$

which can be readily deduced by expanding the Yang-Baxter relation (A.6) to the leading order $\mathcal{O}(\hbar^2)$.

**Sklyanin quadratic Poisson algebra.** In order to deduce the semiclassical limit of the fundamental commutation relation, one has to send $\hbar \to 0$ while keeping $\rho$ and $\hbar S$ constant, thereby fixing the size of classical variables. This is achieved as follows. By expanding Eq. (A.1) to the leading order $\mathcal{O}(\hbar)$ and subsequently replacing commutators with Poisson brackets by virtue of the canonical correspondence principle $^{11}$

$$ [\circ, \circ] \to 2i\hbar\{\circ, \circ\}, \qquad (A.19) $$

the fundamental commutation relation reduces to the Sklyanin quadratic bracket [69]

$$ \{L(\lambda) \otimes L(\lambda')\} = [r(\lambda - \lambda'), L(\lambda) \otimes L(\lambda')], \qquad (A.20) $$

<footnote>$^{10}$There is freedom to choose the overall scale of the $r$-matrix. Here we adopt the common convention [60], 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. [69].</footnote>

<footnote>$^{11}$Normalization of the Poisson bracket is set by the choice of the classical $r$-matrix, cf. Eq. (A.15).</footnote>
with the classical Lax matrix [69]

$$L(\lambda) = 1 S^0 + \frac{1}{i} \sum_a w_a(\lambda) \sigma^a S^a, \qquad (A.21)$$

where we have simultaneously introduced classical ‘Sklyanin variables’ 12

$$\hat{S}^0 \rightarrow \hbar^{-1} S^0, \qquad \hat{S}^a \rightarrow \hbar^{-1} \rho^{-1} S^a, \quad \text{for } a \in \{1, 2, 3\}. \qquad (A.22)$$

Sklyanin quadratic bracket (A.20) provides a compact algebraic representation for the following quadratic Poisson algebra [69]

$$\{S^0, S^3\} = 0, \qquad \{S^0, S^a\} = \rho^{-1} \epsilon_{abc} J_{bc} S^b S^c, \qquad \{S^a, S^b\} = -\rho \epsilon_{abc} S^0 S^c, \qquad (A.23)$$

with classical structure constants $J_{ab} \equiv J_a - J_b$. Poisson algebra (A.23) admits two independent Casimir functions

$$C_0 \equiv \sum_{a=1}^{3} (S^a)^2, \qquad C_1 \equiv (S^0)^2 - (S^3)^2. \qquad (A.24)$$

By prescribing their values $c_0$ and $c_1$, respectively, variables $S^a$ takes values on a two-dimensional non-degenerate Poisson submanifold. By fixing them to

$$c_0 = \sinh^2(\rho), \qquad c_1 = 1, \qquad (A.25)$$

the manifold is topologically equivalent to a two-sphere; $\{S^a\}_{a=1}^3$ lie on a two-sphere, whereas $S^0 > 0$ is uniquely determined by $S^3$ through $C_1$.

Introducing linear combinations $K^\pm \equiv S^0 \pm S^3$ and putting $K \equiv K^+, K^- = C_1$, yielding [114]

$$A_q: \qquad \{K, S^\pm\} = \mp i \rho S^\pm K, \qquad \{S^+, S^-\} = -\frac{i \rho}{2} (K^2 - K^{-2}). \qquad (A.26)$$

The above Poisson relation also follow from quantum algebraic relations (A.11) upon rescaling $\kappa$ as $\kappa = \hbar \rho$, applying the canonical corespondence (A.19), and finally sending $\hbar \rightarrow 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_\ell, y_\ell)$:

(i) $U_1 - U_2$ exchange (35)

$$x_1' = x_1, \qquad (B.1)$$

$$x_2' = x_2 + (u_1 - u_2) \frac{x_1 x_2 v_2 (u_1 + u_2 y_2^2) - x_2^2 (1 + u_1 u_2 v_2^2 y_2^2)}{x_1 u_1 u_2 v_2 (1 - y_2^2) + u_2 x_2 (u_1^2 v_2^2 y_2^2 - 1)}, \qquad (B.2)$$

$$(y_1')^2 = y_1^2 - (u_1^2 - u_2^2) \frac{v_2 x_1 y_1^2 (y_2^2 - 1)}{u_1 x_2 (u_2^2 v_2^2 y_2^2 - 1) + x_1 v_2 u_2^2 (1 - y_2^2)}, \qquad (B.3)$$

$$(y_2')^2 = y_2^2 + (u_1^2 - u_2^2) \frac{v_2 x_1 y_2^2 (y_2^2 - 1)}{u_1 x_2 (u_2^2 v_2^2 y_2^2 - 1) + x_1 v_2 u_1^2 (1 - y_2^2)}, \qquad (B.4)$$

<footnote n="12">¹²Our normalization here slightly differs from the one in Ref. [60] in that $S^a$ for $a \in \{1, 2, 3\}$ were rescaled by a factor of $\rho^{-1}$.</footnote>
and in 'Trotter parametrization' (41)<fnref n="13" />

$$x'_{1} = x_{1}, \qquad (B.5)$$

$$x'_{2} = x_{2} \frac{x_{1}(y_{2}w^{-1} - wy_{2}^{-1}) + x_{2}(v^{-2}y_{2}^{-1} - y_{2}v^{2})}{x_{1}(y_{2} - y_{2}^{-1}) + x_{2}(w^{-1}v^{-2}y_{2}^{-1} - y_{2}wv^{2})}, \qquad (B.6)$$

$$(y'_{1})^{2} = y_{1}^{2} \frac{x_{1}(y_{2} - y_{2}^{-1})w + x_{2}(y_{2}^{-1}v^{-2} - y_{2}v^{2})}{x_{1}(y_{2} - y_{2}^{-1})w^{-1} + x_{2}(y_{2}^{-1}v^{-2} - y_{2}v^{2})}, \qquad (B.7)$$

$$(y'_{2})^{2} = y_{2}^{2} \frac{x_{1}(y_{2} - y_{2}^{-1})w^{-1} + x_{2}(y_{2}^{-1}v^{-2} - y_{2}v^{2})}{x_{1}(y_{2} - y_{2}^{-1})w + x_{2}(y_{2}^{-1}v^{-2} - y_{2}v^{2})}, \qquad (B.8)$$

(ii) $V_1 - V_2$ exchange (37)

$$x'_{1} = x_{1} - (v_{1} - v_{2}) \frac{x_{1}x_{2}u_{1}(v_{2} + v_{1}y_{1}^{2}) - x_{1}^{2}(1 + v_{1}v_{2}u_{1}^{2}y_{1}^{2})}{x_{2}v_{1}v_{2}u_{1}(1 - y_{1}^{2}) + v_{1}x_{1}(u_{1}^{2}v_{2}^{2}y_{1}^{2} - 1)}, \qquad (B.9)$$

$$x'_{2} = x_{2}, \qquad (B.10)$$

$$(y'_{1})^{2} = y_{1}^{2} - (v_{1}^{2} - v_{2}^{2}) \frac{u_{1}x_{2}y_{1}^{2}(y_{1}^{2} - 1)}{v_{2}x_{1}(u_{1}^{2}v_{1}^{2}y_{1}^{2} - 1) + x_{2}u_{1}v_{2}^{2}(1 - y_{1}^{2})}, \qquad (B.11)$$

$$(y'_{2})^{2} = y_{2}^{2} + (v_{1}^{2} - v_{2}^{2}) \frac{u_{1}x_{2}y_{2}^{2}(y_{1}^{2} - 1)}{v_{2}x_{1}(u_{1}^{2}v_{1}^{2}y_{1}^{2} - 1) + x_{2}u_{1}v_{1}^{2}(1 - y_{1}^{2})}, \qquad (B.12)$$

and in 'Trotter parametrization' (41)

$$x'_{1} = x_{1} \frac{x_{1}(y_{1}^{-1}v^{-2} - y_{1}v^{2}) + x_{2}(y_{1}w^{-1} - wy_{1}^{-1})}{x_{1}(y_{1}^{-1}w^{-1}v^{-2} - y_{1}wv^{2}) + x_{2}(y_{1} - y_{1}^{-1})}, \qquad (B.13)$$

$$x'_{2} = x_{2}, \qquad (B.14)$$

$$(y'_{1})^{2} = y_{1}^{2} \frac{x_{1}(y_{1}^{-1}v^{-2} - v^{2}y_{1}) + x_{2}w^{-1}(y_{1} - y_{1}^{-1})}{x_{1}(y_{1}^{-1}v^{-2} - v^{2}y_{1}) + x_{2}w(y_{1} - y_{1}^{-1})}, \qquad (B.15)$$

$$(y'_{2})^{2} = y_{2}^{2} \frac{x_{1}(y_{1}^{-1}v^{-2} - v^{2}y_{1}) + x_{2}w(y_{1} - y_{1}^{-1})}{x_{1}(y_{1}^{-1}v^{-2} - v^{2}y_{1}) + x_{2}w^{-1}(y_{1} - y_{1}^{-1})}, \qquad (B.16)$$

The combined $U-V$ exchange (in Trotter parametrization (41)) gives the two-body map $\Phi$

$$x'_{1} = x_{1} \frac{x_{1}Y_{1}(v, 0) - x_{2}Y_{1}(1, -\tau)}{x_{1}Y_{1}(v, \tau) - x_{2}Y_{1}(1, 0)}, \qquad x'_{2} = x_{2} \frac{x_{1}Y_{2}(1, -\tau) - x_{2}Y_{2}(v, 0)}{x_{1}Y_{2}(1, 0) - x_{2}Y_{2}(v, \tau)}, \qquad (B.17)$$

$$y'_{1}^{2} = y_{1}^{2} \frac{f_{2}^{-}}{f_{1}^{+}}, \qquad y'_{2}^{2} = y_{2}^{2} \frac{f_{1}^{+}}{f_{2}^{-}}, \qquad (B.18)$$

with

$$f_{j}^{\pm} = x_{j \pm 1}^{2} w Y_{j}(1, 0) Y_{j \pm 1}(v, 0) + x_{j}^{2} w^{-1} Y_{j \pm 1}(1, 0) Y_{j}(v, 0) \\ - x_{j \pm 1} x_{j} (Y_{j \pm 1}(v, 0) Y_{j}(v, 0) + Y_{j \pm 1}(1, 0) Y_{j}(1, 0) - y_{j} y_{j \pm 1}^{-1} (w - w^{-1})^{2}), \qquad (B.19)$$

$$Y_{j}(v, \tau) = v^{2}y_{j}w - v^{-2}y_{j}^{-1}w^{-1}. \qquad (B.20)$$

<footnote n="13">¹³Where here and in the following we use $w = e^{i\theta\tau}$.</footnote>
C Isospectral flows and local conservation laws

A key implication of the discrete zero-curvature condition (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_l}; λ)$, defined as a right-to-left ordered product of Lax matrices along a zig-zag path (i.e. initial sawtooth) of length L

$$T_{\tau}(\lambda) = L(S_{L}; \lambda^{+})L(S_{L-1}; \lambda^{-})\cdots L(S_{2}; \lambda^{+})L(S_{1}; \lambda^{-}), \qquad (C.1)$$

where $λ^± ≡ λ±τ/2$. Here and subsequently we leave the dependence on anisotropy $\rho$ implicit. By tracing out the auxiliary $\mathbb{C}^2$ space, we obtain the transfer function

$$t_{\tau}(\lambda) = \mathrm{Tr} T_{\tau}(\lambda). \qquad (C.2)$$

**Isospectral evolution.** By virtue of the zero-curvature condition (3) and cyclic property of the trace, it follows that

$$t_{\tau}(\lambda) \circ \Phi_{\tau}^{\text{even}} = t_{-\tau}(\lambda), \qquad t_{-\tau}(\lambda) \circ \Phi_{\tau}^{\text{odd}} = t_{\tau}(\lambda), \qquad (C.3)$$

which in combination imply invariance under the time propagation (6)

$$t_{\tau}(\lambda) \circ \Phi_{\tau}^{\text{full}} = t_{\tau}(\lambda). \qquad (C.4)$$

The quadratic bracket (A.20) can be immediately lifted onto the level of monodromies

$$\{T_{\tau}(\lambda)^{\otimes}, T_{\tau}(\mu)\} = [r(\lambda - \mu), T_{\tau}(\lambda) \otimes T_{\tau}(\mu)]. \qquad (C.5)$$

By subsequently taking the trace over $\mathbb{C}^2$, one readily finds that transfer functions Poisson commute,

$$\{t_{\tau}(\lambda), t_{\tau}(\lambda')\} = 0, \qquad (C.6)$$

for any pair of complex spectral parameters $λ, λ' \in \mathbb{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 \to \infty$, 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

$$\overline{L(\lambda)} = \sigma^2 L(\overline{\lambda}) \sigma^2, \qquad L^T(\lambda) = \sigma^2 L(-\lambda) \sigma^2, \qquad (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 separate infinite towers of local conserved quantities, denoting them by $Q_{k;\tau}^{\pm}$, using the standard approach (see, for example, Ref. [60]) by expanding the logarithmic derivatives of the modulus of $t_τ(λ)$ around either of the two distinguished points $λ_0^{\pm} = i \mp \tau/2$ at which the Lax matrix degenerates (becoming a projector of rank one),

$$\mathrm{Det} L(\lambda_0^{\pm}) = 0, \qquad L(\lambda_0^{\pm}) = |\alpha_{\pm}\rangle \langle \beta_{\pm}|, \qquad (C.8)$$
yielding two independent sequences of local charges

$$Q_{k;\tau}^{\pm} = \left. \frac{\partial^k}{\partial \lambda^k} \log |t_{\tau}(\lambda)|^2 \right|_{\lambda=\lambda_0^{\pm}}. \qquad (C.9)$$

Taking the continuous time limit $\tau \to 0$, both branches merge into the local conservation laws of the anisotropic lattice Landau-Lifshitz model.

**Parity symmetry.** The transfer function $t_{\tau}(\lambda)$ is $\mathbb{Z}_2$-invariant under the spin-reversal transformation ($\pi$-rotation of the canonical spin $\mathbf{S}$ around the $x$-axis)

$$\mathcal{F}_1: \quad S^3 \mapsto -S^3, \quad S^{\pm} \mapsto S^{\mp}, \qquad (C.10)$$

implying $\mathcal{F}_1(\mathcal{K}) = \mathcal{K}^{-1}$. At the level of the Lax matrix, $\mathcal{F}_1$ represents conjugation by Pauli matrix $\sigma^1$,

$$\mathcal{F}_1: \quad L(\lambda) \mapsto \sigma^1 L(\lambda) \sigma^1, \quad T(\lambda) \mapsto \sigma^1 T(\lambda) \sigma^1, \qquad (C.11)$$

which implies invariance $\mathcal{F}_1(t_{\tau}(\lambda)) = t_{\tau}(\lambda)$ and hence charges $Q_{k;\tau}^{\pm}$ are all invariant under $\mathcal{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 $\rho$, discrete time-step $\tau$, and lattice spacing $\Delta$ (which we have kept fixed $\Delta = 1$ throughout the text). The full set of models which descend from $\Phi_{\tau,\rho}$ can be therefore arranged on the vertices of a cube, as depicted in Figure 10. We briefly describe them below.<fnref n="14" />

<figure><img src="image_13.png" /><figcaption>Figure 10: Reductions of the anisotropic discrete-time propagator $\Phi_{\tau,\rho}$ (upper-left corner), obtained via taking various combinations of $\rho, \Delta, \tau \to 0$ limits. The complete reduction obtained by taking all three limit yields the isotropic Landau-Lifshitz field theory $\mathcal{H}_{LL}$ (lower-right corner).</figcaption></figure>

<footnote n="14">¹⁴We shall omit discussing $\Psi_{\tau}^{\rho}$ and $\Psi_{\tau}$, which correspond to (anisotropic) 'kicked' field theories in discrete time.</footnote>
**Φ<sub>τ</sub> - isotropic discrete space-time model.** By letting the anisotropy parameter to zero, $\rho \to 0$, the anisotropy Lax matrix given by Eq. (15) becomes the Lax matrix of the isotropic lattice Landau-Lifshitz model [69, 115]

$$L^{(0)}(\lambda, \mathbf{S}) = 1 + \frac{1}{i\lambda} \mathbf{S} \cdot \boldsymbol{\sigma}, \qquad (D.1)$$

where $\boldsymbol{\sigma} = (\sigma^1, \sigma^2, \sigma^3)$ is a vector of Pauli matrices. The corresponding discrete zero-curvature condition around a plaquette,

$$L_2^{(0)}(\lambda + \tau/2)L_1^{(0)}(\lambda - \tau/2) = L_2^{(0)\prime}(\lambda - \tau/2)L_1^{(0)}(\lambda + \tau/2), \qquad (D.2)$$

generates a two-body symplectic propagator [55]

$$\Phi_{\tau}(\mathbf{S}_1, \mathbf{S}_2) = \frac{1}{s^2 + \tau^2} (s^2 \mathbf{S}_1 + \tau^2 \mathbf{S}_2 + \tau \mathbf{S}_1 \times \mathbf{S}_2, s^2 \mathbf{S}_2 + \tau^2 \mathbf{S}_1 + \tau \mathbf{S}_2 \times \mathbf{S}_1), \qquad (D.3)$$

$$s^2 = \frac{1}{2}(1 + \mathbf{S}_1 \cdot \mathbf{S}_2). \qquad (D.4)$$

**$H_{LLL}^{\rho}$ - uniaxially anisotropic lattice Landau-Lifshitz model.** By sending $\tau \to 0$, the discrete-time propagator becomes the equation of motion in continuous time,

$$\frac{d}{dt} S_{\ell}^{a} = \{S^{a}, H_{LLL}^{\rho}\}, \qquad (D.5)$$

generated by a Hamiltonian

$$H_{LLL}^{\rho} \approx \sum_{\ell=1}^{L} \log \left[ \sinh^2 \rho S_{\ell}^{0} S_{\ell+1}^{0} + S_{\ell}^{1} S_{\ell+1}^{1} + S_{\ell}^{2} S_{\ell+1}^{2} + \cosh^2 \rho S_{\ell}^{3} S_{\ell+1}^{3} \right]. \qquad (D.6)$$

**$H_{LLL}$ - isotropic lattice Landau-Lifshitz model.** Starting either from Eq. (D.4) and sending $\tau \to 0$, or from Eq. (D.5) and sending $\rho \to 0$, one ends up with the Hamiltonian dynamics of the isotropic lattice Landau-Lifshitz model. In terms of canonical spins $\mathbf{S}$, the equations of motion read [115]

$$\frac{d}{dt} \mathbf{S}_{\ell} = \{\mathbf{S}_{\ell}, H_{LLL}\} = \frac{\mathbf{S}_{\ell} \times \mathbf{S}_{\ell-1}}{1 + \mathbf{S}_{\ell} \cdot \mathbf{S}_{\ell-1}} + \frac{\mathbf{S}_{\ell} \times \mathbf{S}_{\ell+1}}{1 + \mathbf{S}_{\ell} \cdot \mathbf{S}_{\ell+1}}, \qquad (D.7)$$

with a Hamiltonian of the form $H_{LLL} = \sum_{\ell=1}^{L} \log(1 + \mathbf{S}_{\ell} \cdot \mathbf{S}_{\ell+1})$.

**$H_{LLL}^{\rho}$ - anisotropic Landau-Lifshitz field theory.** To obtain the full continuum limit, we reinstate the lattice spacing $\Delta$ and subsequently retain smooth lattice spin configurations $\mathbf{S}_{\ell} = \mathbf{S}(x = \Delta\ell)$ by keeping only variations at the leading non-trivial order $O(\Delta^2)$:

$$\mathbf{S}_{\ell \pm 1} = \mathbf{S}(x) \pm \Delta \partial_x \mathbf{S} + \frac{1}{2} \Delta^2 \partial_x^2 \mathbf{S} + \mathcal{O}(\Delta^3). \qquad (D.8)$$

Upon taking the continuum limit $\Delta \to 0$, both the Sklyanin variables and the anisotropy parameter have to be accordingly rescaled by the lattice spacing, namely

$$S^0 = 1 + \mathcal{O}(\Delta^2), \qquad S^j = \Delta S^a + \mathcal{O}(\Delta^2) \quad \text{for } a = 1, 2, 3, \qquad (D.9)$$

and $\rho \to \Delta\rho$. This way, one recovers the anisotropic Landau-Lifshitz field theory [59, 60]

$$H_{LL}^{\rho} = \frac{1}{2} \int_{0}^{L} dx \left[ (\partial_{x} \mathbf{S})^{2} - \mathbf{S} \cdot \mathcal{J} \mathbf{S} \right], \qquad \mathcal{J} = \text{diag}(J_{1}, J_{2}, J_{3}), \qquad (D.10)$$

(where $J_1 = J_2$ and $J_3 - J_1 = \rho^2$) with the equation of motion [59, 60]

$$\partial_t \mathbf{S}(x) = \{\mathbf{S}(x), H_{LL}^{\rho}\} = \mathbf{S} \times \partial_x^2 \mathbf{S} + \mathbf{S} \times \mathcal{J} \mathbf{S}. \qquad (D.11)$$
$\mathcal{H}_{LL}$ – isotropic Landau-Lifshitz field theory. By finally sending $\rho \to 0$, Eq. (D.10) becomes the isotropic Landau-Lifshitz field theory [58, 60]

$$ \mathcal{H}_{LL} = \frac{1}{2} \int_{0}^{L} dx (\partial_{x} \mathbf{S})^{2}, \quad \partial_{t} \mathbf{S} = \{\mathbf{S}, \mathcal{H}_{LL}\} = \mathbf{S} \times \partial_{x}^{2} \mathbf{S}. \qquad (D.12) $$

References

[1] G. Vidal, *Efficient classical simulation of slightly entangled quantum computations*, Phys. Rev. Lett. **91**, 147902 (2003), doi:10.1103/PhysRevLett.91.147902.

[2] F. Verstraete, J. J. García-Ripoll and J. I. Cirac, *Matrix product density operators: Simulation of finite-temperature and dissipative systems*, Phys. Rev. Lett. **93**, 207204 (2004), doi:10.1103/PhysRevLett.93.207204.

[3] D. A. Abanin, E. Altman, I. Bloch and M. Serbyn, *Colloquium: Many-body localization, thermalization, and entanglement*, Rev. Mod. Phys. **91**, 021001 (2019), doi:10.1103/RevModPhys.91.021001.

[4] R. Nandkishore and D. A. Huse, *Many-body localization and thermalization in quantum statistical mechanics*, Annu. Rev. Condens. Matter Phys. **6**, 15 (2015), doi:10.1146/annurev-conmatphys-031214-014726.

[5] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg and M. Žnidarič, *Finite-temperature transport in one-dimensional quantum lattice models*, Rev. Mod. Phys. **93**, 025003 (2021), doi:10.1103/RevModPhys.93.025003.

[6] L. D’Alessio, Y. Kafri, A. Polkovnikov and M. Rigol, *From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics*, Advances in Physics **65**, 239 (2016), doi:10.1080/00018732.2016.1198134.

[7] I. Bloch, J. Dalibard and W. Zwerger, *Many-body physics with ultracold gases*, Rev. Mod. Phys. **80**, 885 (2008), doi:10.1103/revmodphys.80.885.

[8] I. Bloch, J. Dalibard and S. Nascimbène, *Quantum simulations with ultracold quantum gases*, Nat. Phys. **8**, 267 (2012), doi:10.1038/nphys2259.

[9] T. Langen et al., *Experimental observation of a generalized Gibbs ensemble*, Science **348**, 207 (2015), doi:10.1126/science.1257026.

[10] M. Schemmer, I. Bouchoule, B. Doyon and J. Dubail, *Generalized hydrodynamics on an atom chip*, Phys. Rev. Lett. **122**, 090601 (2019), doi:10.1103/physrevlett.122.090601.

[11] P. N. Jepsen, J. Amato-Grill, I. Dimitrova, W. Wei Ho, E. Demler and W. Ketterle, *Spin transport in a tunable Heisenberg model realized with ultracold atoms*, Nature **588**, 403 (2020), doi:10.1038/s41586-020-3033-y.

[12] Y. Le, N. Malvania, Y. Zhang, J. Dubail, M. Rigol and D. Weiss, *Generalized hydrodynamics in strongly interacting 1D Bose gases (I)*, Science **373**, 1129 (2021), doi:10.1126/science.abf0147.

[13] B. Doyon, J. Dubail, R. Konik and T. Yoshimura, *Large-scale description of interacting one-dimensional Bose gases: Generalized hydrodynamics supersedes conventional hydrodynamics*, Phys. Rev. Lett. **119**, 195301 (2017), doi:10.1103/physrevlett.119.195301.
[14] G. Misguich, K. Mallick and P. L. Krapivsky, *Dynamics of the spin-1/2 Heisenberg chain initialized in a domain-wall state*, Phys. Rev. B **96**, 195151 (2017), doi:10.1103/PhysRevB.96.195151.

[15] M. Ljubotina, M. Žnidarič and T. Prosen, *A class of states supporting diffusive spin dynamics in the isotropic Heisenberg model*, J. Phys. A: Math. Theor. **50**, 475002 (2017), doi:10.1088/1751-8121/aa8bdc.

[16] B. Doyon, T. Yoshimura and J.-S. Caux, *Soliton gases and generalized hydrodynamics*, Phys. Rev. Lett. **120**, 045301 (2018), doi:10.1103/physrevlett.120.045301.

[17] B. Doyon and T. Yoshimura, *A note on generalized hydrodynamics: inhomogeneous fields and other concepts*, SciPost Phys. **2**, 014 (2017), doi:10.21468/SciPostPhys.2.2.014.

[18] A. Bastianello, J. De Nardis and A. De Luca, *Generalized hydrodynamics with dephasing noise*, Phys. Rev. B **102**, 161110 (2020), doi:10.1103/physrevb.102.161110.

[19] J.-S. Caux, B. Doyon, J. Dubail, R. Konik and T. Yoshimura, *Hydrodynamics of the interacting Bose gas in the Quantum Newton Cradle setup*, SciPost Phys. **6**, 070 (2019), doi:10.21468/scipostphys.6.6.070.

[20] B. Bertini, L. Piroli and M. Kormos, *Transport in the sine-Gordon field theory: From generalized hydrodynamics to semiclassical*, Phys. Rev. B **100**, 035108 (2019), doi:10.1103/physrevb.100.035108.

[21] A. Urichuk, Y. Oez, A. Klümper and J. Sirker, *The spin Drude weight of the XXZ chain and generalized hydrodynamics*, SciPost Phys. **6**, 005 (2019), doi:10.21468/SciPostPhys.6.1.005.

[22] M. Borsi, B. Pozsgay and L. Pristyák, *Current operators in Bethe ansatz and generalized hydrodynamics: An exact quantum-classical correspondence*, Phys. Rev. X **10**, 011054 (2020), doi:10.1103/physrevx.10.011054.

[23] A. C. Cassidy, C. W. Clark and M. Rigol, *Generalized thermalization in an integrable lattice system*, Phys. Rev. Lett. **106**, 140405 (2011), doi:10.1103/physrevlett.106.140405.

[24] L. Vidmar and M. Rigol, *Generalized Gibbs ensemble in integrable lattice models*, J. Stat. Mech. 064007 (2016), doi:10.1088/1742-5468/2016/06/064007.

[25] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler and T. Prosen, *Complete generalized Gibbs ensembles in an interacting theory*, Phys. Rev. Lett. **115**, 157201 (2015), doi:10.1103/PhysRevLett.115.157201.

[26] F. H. L. Essler and M. Fagotti, *Quench dynamics and relaxation in isolated integrable quantum spin chains*, J. Stat. Mech. 064002 (2016), doi:10.1088/1742-5468/2016/06/064002.

[27] B. Bertini, M. Collura, J. De Nardis and M. Fagotti, *Transport in out-of-equilibrium XXZ chains: Exact profiles of charges and currents*, Phys. Rev. Lett. **117**, 207201 (2016), doi:10.1103/PhysRevLett.117.207201.

[28] O. A. Castro-Alvaredo, B. Doyon and T. Yoshimura, *Emergent hydrodynamics in integrable quantum systems out of equilibrium*, Phys. Rev. X **6**, 041065 (2016), doi:10.1103/PhysRevX.6.041065.
[29] V. B. Bulchandani, R. Vasseur, C. Karrasch and J. E. Moore, *Solvable hydro-dynamics of quantum integrable systems*, Phys. Rev. Lett. **119**, 220604 (2017), doi:10.1103/physrevlett.119.220604.

[30] M. Ljubotina, M. Žnidarič and T. Prosen, *Spin diffusion from an inhomogeneous quench in an integrable system*, Nat. Commun. **8**, 16117 (2017), doi:10.1038/ncomms16117.

[31] V. B. Bulchandani, R. Vasseur, C. Karrasch and J. E. Moore, *Bethe-Boltzmann hydrodynamics and spin transport in the XXZ chain*, Phys. Rev. B **97**, 045407 (2018), doi:10.1103/physrevb.97.045407.

[32] M. Ljubotina, M. Žnidarič and T. Prosen, *Kardar-Parisi-Zhang physics in the quantum Heisenberg magnet*, Phys. Rev. Lett. **122**, 210602 (2019), doi:10.1103/physrevlett.122.210602.

[33] M. Dupont and J. E. Moore, *Universal spin dynamics in infinite-temperature one-dimensional quantum magnets*, Phys. Rev. B **101**, 121106 (2020), doi:10.1103/physrevb.101.121106.

[34] S. Gopalakrishnan and R. Vasseur, *Kinetic theory of spin diffusion and superdiffusion in XXZ spin chains*, Phys. Rev. Lett. **122**, 127202 (2019), doi:10.1103/physrevlett.122.127202.

[35] S. Gopalakrishnan, R. Vasseur and B. Ware, *Anomalous relaxation and the high-temperature structure factor of XXZ spin chains*, Proc. Natl. Acad. Sci. USA **116**, 16250 (2019), doi:10.1073/pnas.1906914116.

[36] J. De Nardis, M. Medenjak, C. Karrasch and E. Ilievski, *Universality classes of spin transport in one-dimensional isotropic magnets: The onset of logarithmic anomalies*, Phys. Rev. Lett. **124**, 210605 (2020), doi:10.1103/physrevlett.124.210605.

[37] U. Agrawal, S. Gopalakrishnan, R. Vasseur and B. Ware, *Anomalous low-frequency conductivity in easy-plane XXZ spin chains*, Phys. Rev. B **101**, 224415 (2020), doi:10.1103/PhysRevB.101.224415.

[38] V. Alba, B. Bertini, M. Fagotti, L. Piroli and P. Ruggiero, *Generalized-hydrodynamic approach to inhomogeneous quenches: Correlations, entanglement and quantum effects* (2021), arXiv:2104.00656.

[39] A. Cortés Cubero, T. Yoshimura and H. Spohn, *Form factors and generalized hydrodynamics for integrable systems* (2021), arXiv:2104.04951.

[40] J. De Nardis, B. Doyon, M. Medenjak and M. Panfil, *Correlation functions and transport coefficients in generalised hydrodynamics* (2021), arXiv:2104.04462.

[41] V. B. Bulchandani, S. Gopalakrishnan and E. Ilievski, *Superdiffusion in spin chains*, J. Stat. Mech. 084001 (2021), doi:10.1088/1742-5468/ac12c7.

[42] B. Buča, K. Klobas and T. Prosen, *Rule 54: exactly solvable model of nonequilibrium statistical mechanics*, J. Stat. Mech. 074001 (2021), doi:10.1088/1742-5468/ac096b.

[43] M. Borsi, B. Pozsgay and L. Prstyák, *Current operators in integrable models: A review* (2021), arXiv:2103.12160.

[44] A. Bastianello, B. Doyon, G. Watts and T. Yoshimura, *Generalized hydrodynamics of classical integrable field theory: the sinh-Gordon model*, SciPost Phys. **4**, 045 (2018), doi:10.21468/scipostphys.4.6.045.
[45] B. Doyon, *Generalized hydrodynamics of the classical Toda system*, J. Math. Phys. **60**, 073302 (2019), doi:10.1063/1.5096892.

[46] O. Gamayun, Y. Miao and E. Ilievski, *Domain-wall dynamics in the Landau-Lifshitz magnet and the classical-quantum correspondence for spin transport*, Phys. Rev. B **99**, 140301 (2019), doi:10.1103/physrevb.99.140301.

[47] H. Spohn, *Collision rate ansatz for the classical Toda lattice*, Phys. Rev. E **101**, 060103 (2020), doi:10.1103/physreve.101.060103.

[48] H. Spohn, *Ballistic space-time correlators of the classical toda lattice*, J. Phys. A: Math. Theor. **53**, 265004 (2020), doi:10.1088/1751-8121/ab91d5.

[49] A. Kuniba, G. Misguich and V. Pasquier, *Generalized hydrodynamics in box-ball system*, J. Phys. A: Math. Theor. **53**, 404001 (2020), doi:10.1088/1751-8121/abadb9.

[50] D. A. Croydon and M. Sasada, *Generalized hydrodynamic limit for the box-ball system*, Commun. Math. Phys. **383**, 427 (2020), doi:10.1007/s00220-020-03914-x.

[51] A. Kuniba, G. Misguich and V. Pasquier, *Generalized hydrodynamics in complete box-ball system for $U_q(\hat{s}\hat{l}_n)$*, SciPost Phys. **10**, 095 (2021), doi:10.21468/scipostPhys.10.4.095.

[52] M. Vanicat, L. Zadnik and T. Prosen, *Integrable trotterization: Local conservation laws and boundary driving*, Phys. Rev. Lett. **121**, 030606 (2018), doi:10.1103/PhysRevLett.121.030606.

[53] Ž. Krajnik, E. Ilievski and T. Prosen, *Integrable matrix models in discrete space-time*, SciPost Phys. **9**, 038 (2020), doi:10.21468/scipostphys.9.3.038.

[54] K. Klobas, M. Medenjak and T. Prosen, *Exactly solvable deterministic lattice model of crossover between ballistic and diffusive transport*, J. Stat. Mech. 123202 (2018), doi:10.1088/1742-5468/aae853.

[55] Ž. Krajnik and T. Prosen, *Kardar-Parisi-Zhang physics in integrable rotationally symmetric dynamics on discrete space-time lattice*, J. Stat. Phys. **179**, 110 (2020), doi:10.1007/s10955-020-02523-1.

[56] V. Gritsev and A. Polkovnikov, *Integrable floquet dynamics*, SciPost Phys. **2**, 021 (2017), doi:10.21468/scipostphys.2.3.021.

[57] M. Ljubotina, L. Zadnik and T. Prosen, *Ballistic spin transport in a periodically driven integrable quantum system*, Phys. Rev. Lett. **122**, 150605 (2019), doi:10.1103/physrevlett.122.150605.

[58] L. A. Takhtajan, *Integration of the continuous Heisenberg spin chain through the inverse scattering method*, Phys. Lett. A **64**, 235 (1977), doi:10.1016/0375-9601(77)90727-7.

[59] E. Sklyanin, *On complete integrability of the Landau-Lifshitz equation* (1979), LOMI-79-3.

[60] L. D. Faddeev and L. A. Takhtajan, *Hamiltonian methods in the theory of solitons*, Springer Berlin Heidelberg, ISBN 9783540698432 (1987), doi:10.1007/978-3-540-69969-9.

[61] R. Hirota, *Bilinearization of soliton equations*, J. Phys. Soc. Jpn. **51**, 323 (1982), doi:10.1143/jpsj.51.323.
[62] P. D. Lax, *Integrals of nonlinear equations of evolution and solitary waves*, Comm. Pure Appl. Math. **21**, 467 (1968), doi:10.1002/cpa.3160210503.

[63] M. J. Ablowitz, D. J. Kaup, A. C. Newell and H. Segur, *The inverse scattering transform-Fourier analysis for nonlinear problems*, Stud. Appl. Math. **53**, 249 (1974), doi:10.1002/sapm1974534249.

[64] M. J. Ablowitz and H. Segur, *Solitons and the inverse scattering transform*, SIAM, ISBN 9780898711745 (1981), doi:10.1137/1.9781611970883.

[65] O. Babelon, D. Bernard and M. Talon, *Introduction to classical integrable systems*, Cambridge University Press, Cambridge, ISBN 9780521822671 (2003), doi:10.1017/cbo9780511535024.

[66] J. Hietarinta, N. Joshi and F. W. Nijhoff, *Discrete systems and integrability*, Cambridge University Press, Cambridge, ISBN 9781107337411 (2016), doi:10.1017/cbo9781107337411.

[67] J. Moser and A. P. Veselov, *Discrete versions of some classical integrable systems and factorization of matrix polynomials*, Commun. Math. Phys. **139**, 217 (1991), doi:10.1007/bf02352494.

[68] A. P. Veselov, *Yang-Baxter maps and integrable dynamics*, Phys. Lett. A **314**, 214 (2003), doi:10.1016/s0375-9601(03)00915-0.

[69] E. K. Sklyanin, *Some algebraic structures connected with the Yang-Baxter equation*, Funct. Anal. Appl. **16**, 263 (1982), doi:10.1007/bf01077848.

[70] T. Prosen and B. Žunkovič, *Macroscopic diffusive transport in a microscopically integrable Hamiltonian system*, Phys. Rev. Lett. **111**, 040602 (2013), doi:10.1103/PhysRevLett.111.040602.

[71] A. Das, M. Kulkarni, H. Spohn and A. Dhar, *Kardar-Parisi-Zhang scaling for an integrable lattice Landau-Lifshitz spin chain*, Phys. Rev. E **100**, 042116 (2019), doi:10.1103/PhysRevE.100.042116.

[72] E. K. Sklyanin, *Poisson structure of a periodic classical XYZ-chain*, J. Math. Sci. **46**, 1664 (1989), doi:10.1007/bf01099198.

[73] L. D. Faddeev, N. Yu. Reshetikhin and L. A. Takhtajan, *Quantization of Lie Groups and Lie Algebras*, in Algebraic Analysis, Elsevier, ISBN 9780124004658 (1988), doi:10.1016/b978-0-12-400465-8.50019-5.

[74] V. G. Drinfel'd, *Quantum groups*, J. Math. Sci. **41**, 898 (1988), doi:10.1007/bf01247086.

[75] S. M. Khoroshkin and V. N. Tolstoy, *Yangian double*, Lett. Math. Phys. **36**, 373 (1996), doi:10.1007/bf00714404.

[76] S. Khoroshkin and Z. Tsuboi, *The universal R-matrix and factorization of the L-operators related to the Baxter Q-operators*, J. Phys. A: Math. Theor. **47**, 192003 (2014), doi:10.1088/1751-8113/47/19/192003.

[77] R. J. Baxter, *Exactly solved models in statistical mechanics*, Elsevier (2016).

[78] S. E. Derkachov, *Factorization of the R-matrix. I*, J. Math. Sci. **143**, 2773 (2007), doi:10.1007/s10958-007-0164-8.
[79] S. E. Derkachov, *Factorization of the R-matrix. II*, J. Math. Sci. **143**, 2791 (2007), doi:10.1007/s10958-007-0165-7.

[80] S. Derkachov, D. Karakhanyan and R. Kirschner, *Baxter Q-operators of the XXZ chain and R-matrix factorization*, Nucl. Phys. B **738**, 368 (2006), doi:10.1016/j.nuclphysb.2005.12.015.

[81] V. V. Bazhanov, T. Łukowski, C. Meneghelli and M. Staudacher, *A shortcut to the Q-operator*, J. Stat. Mech. P11002 (2010), doi:10.1088/1742-5468/2010/11/p11002.

[82] R. Frassek, T. Łukowski, C. Meneghelli and M. Staudacher, *Oscillator construction of Q-operators*, Nucl. Phys. B **850**, 175 (2011), doi:10.1016/j.nuclphysb.2011.04.008.

[83] O. Babelon, K. K. Kozlowski and V. Pasquier, *Baxter operator and Baxter equation for q-Toda and Toda2 chains*, Rev. Math. Phys. **30**, 1840003 (2018), doi:10.1142/s0129055x18400032.

[84] A. P. Veselov, *Yang–Baxter maps and integrable dynamics*, Phys. Lett. A **314**, 214 (2003), doi:10.1016/S0375-9601(03)00915-0.

[85] V. V. Bazhanov and S. M. Sergeev, *Yang–Baxter maps, discrete integrable equations and quantum groups*, Nucl. Phys. B **926**, 509 (2018), doi:10.1016/j.nuclphysb.2017.11.017.

[86] Z. Tsuboi, *Quantum groups, Yang–Baxter maps and quasi-determinants*, Nucl. Phys. B **926**, 200 (2018), doi:10.1016/j.nuclphysb.2017.11.005.

[87] B. Doyon, *Lecture notes on generalised hydrodynamics*, SciPost Phys. Lect. Notes 18 (2020), doi:10.21468/scipostphyslectnotes.18.

[88] M. Žnidarič, *Spin transport in a one-dimensional anisotropic Heisenberg model*, Phys. Rev. Lett. **106**, 220601 (2011), doi:10.1103/PhysRevLett.106.220601.

[89] E. Ilievski, J. De Nardis, M. Medenjak and T. Prosen, *Superdiffusion in one-dimensional quantum lattice models*, Phys. Rev. Lett. **121**, 230602 (2018), doi:10.1103/PhysRevLett.121.230602.

[90] J. De Nardis, M. Medenjak, C. Karrasch and E. Ilievski, *Anomalous spin diffusion in one-dimensional antiferromagnets*, Phys. Rev. Lett. **123**, 186601 (2019), doi:10.1103/PhysRevLett.123.186601.

[91] V. B. Bulchandani, *Kardar–Parisi–Zhang universality from soft gauge modes*, Phys. Rev. B **101**, 041411 (2020), doi:10.1103/physrevb.101.041411.

[92] J. De Nardis, S. Gopalakrishnan, E. Ilievski and R. Vasseur, *Superdiffusion from emergent classical solitons in quantum spin chains*, Phys. Rev. Lett. **125**, 070601 (2020), doi:10.1103/PhysRevLett.125.070601.

[93] E. Ilievski, J. De Nardis, S. Gopalakrishnan, R. Vasseur and B. Ware, *Superuniversality of superdiffusion*, Phys. Rev. X **11**, 031023 (2021), doi:10.1103/PhysRevX.11.031023.

[94] M. Kardar, G. Parisi and Y.-C. Zhang, *Dynamic scaling of growing interfaces*, Phys. Rev. Lett. **56**, 889 (1986), doi:10.1103/PhysRevLett.56.889.

[95] E. Ilievski and J. De Nardis, *Microscopic origin of ideal conductivity in integrable quantum models*, Phys. Rev. Lett. **119**, 020602 (2017), doi:10.1103/PhysRevLett.119.020602.
[96] P. Mazur, *Non-ergodicity of phase functions in certain systems*, Physica **43**, 533 (1969), doi:10.1016/0031-8914(69)90185-2.

[97] M. Suzuki, *Ergodicity, constants of motion, and bounds for susceptibilities*, Physica **51**, 277 (1971), doi:10.1016/0031-8914(71)90226-6.

[98] E. Ilievski and T. Prosen, *Thermodynamic bounds on Drude weights in terms of almost-conserved quantities*, Commun. Math. Phys. **318**, 809 (2013), doi:10.1007/s00220-012-1599-4.

[99] B. Doyon, *Diffusion and superdiffusion from hydrodynamic projection* (2019), arXiv:1912.01551.

[100] E. Ilievski, M. Medenjak and T. Prosen, *Quasilocal conserved operators in the isotropic Heisenberg spin-1/2 chain*, Phys. Rev. Lett. **115**, 120601 (2015), doi:10.1103/physrevlett.115.120601.

[101] E. Ilievski, M. Medenjak, T. Prosen and L. Zadnik, *Quasilocal charges in integrable lattice systems*, J. Stat. Mech. 064008 (2016), doi:10.1088/1742-5468/2016/06/064008.

[102] T. Prosen, *Open XXZ spin chain: Nonequilibrium steady state and a strict bound on ballistic transport*, Phys. Rev. Lett. **106**, 217206 (2011), doi:10.1103/PhysRevLett.106.217206.

[103] T. Prosen and E. Ilievski, *Families of quasilocal conservation laws and quantum spin transport*, Phys. Rev. Lett. **111**, 057203 (2013), doi:10.1103/PhysRevLett.111.057203.

[104] T. Prosen, *Quasilocal conservation laws in XXZ spin-1/2 chains: Open, periodic and twisted boundary conditions*, Nucl. Phys. B **886**, 1177 (2014), doi:10.1016/j.nuclphysb.2014.07.024.

[105] R. G. Pereira, V. Pasquier, J. Sirker and I. Affleck, *Exactly conserved quasilocal operators for the XXZ spin chain*, J. Stat. Mech. P09037 (2014), doi:10.1088/1742-5468/2014/09/p09037.

[106] B. Doyon and H. Spohn, *Drude weight for the Lieb-Liniger Bose gas*, SciPost Phys. **3**, 039 (2017), doi:10.21468/scipostphys.3.6.039.

[107] E. Ilievski and J. De Nardis, *Ballistic transport in the one-dimensional Hubbard model: The hydrodynamic approach*, Phys. Rev. B **96**, 081118 (2017), doi:10.1103/physrevb.96.081118.

[108] J. De Nardis, D. Bernard and B. Doyon, *Hydrodynamic diffusion in integrable systems*, Phys. Rev. Lett. **121**, 160603 (2018), doi:10.1103/physrevlett.121.160603.

[109] S. Gopalakrishnan, D. A. Huse, V. Khemani and R. Vasseur, *Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems*, Phys. Rev. B **98**, 220303 (2018), doi:10.1103/physrevb.98.220303.

[110] J. De Nardis, D. Bernard and B. Doyon, *Diffusion in generalized hydrodynamics and quasiparticle scattering*, SciPost Phys. **6**, 049 (2019), doi:10.21468/scipostphys.6.4.049.

[111] V. E. Zakharov, *The inverse scattering method*, in Topics in Current Physics, Springer Berlin Heidelberg, ISBN 978-3-642-81448-8 (1980), doi:10.1007/978-3-642-81448-8_7.
[112] V. E. Korepin, N. M. Bogoliubov and A. G. Izergin, *Quantum inverse scattering method and correlation functions*, Cambridge University Press, ISBN 9780521373203 (1993), doi:10.1017/CBO9780511628832.

[113] E. K. Sklyanin, *Classical limits of the SU(2)-invariant solutions of the Yang-Baxter equation*, J. Math. Sci. **40**, 93 (1988), doi:10.1007/bf01084941.

[114] E. K. Sklyanin, *Poisson structure of classical XXZ-chain*, J. Math. Sci. **46**, 2104 (1989), doi:10.1007/bf01096094.

[115] Y. Ishimori, *An integrable classical spin chain*, J. Phys. Soc. Jpn. **51**, 3417 (1982), doi:10.1143/jpsj.51.3417.