Integrable Matrix Models in Discrete Space-Time

We introduce a class of integrable dynamical systems of interacting classical matrix-valued fields propagating on a discrete space-time lattice, realized as many-body circuits built from elementary symplectic two-body maps. The models provide an efficient integrable Trotterization of non-relativistic $\sigma$-models with complex Grassmannian manifolds as target spaces, including, as special cases, the higher-rank analogues of the Landau--Lifshitz field theory on complex projective spaces. As an application, we study transport of Noether charges in canonical local equilibrium states. We find a clear signature of superdiffusive behavior in the Kardar--Parisi--Zhang universality class, irrespectively of the chosen underlying global unitary symmetry group and the quotient structure of the compact phase space, providing a strong indication of superuniversal physics.


Introduction
Explaining how macroscopic laws of matter emerge from microscopic reversible dynamics is one the central problems of modern theoretical physics which is still quite far from being settled. One source of difficulties is that dynamical systems which comprise many interacting degrees of freedom are rarely amenable to exact analytic treatment and explicit closed-form solutions are an exception. To make matter worse, numerical simulations of thermodynamic systems out of equilibrium becomes quickly inaccessible at large times owing to exponential growth of required resources, even in one spatial dimension where state-of-the-art methods based on matrix-product states are available. Integrable models provide an opportunity to mitigate some of these issues by providing an ideal theoretical playground and address some key question of statistical physics with a high level of rigour. In spite of a long-lasting progress in the field of classical [1][2][3][4][5][6] and quantum integrability [7][8][9][10][11][12][13][14][15][16], the ultimate hope to obtain explicit solutions to various nonequilibrium problems has not materialized yet, and even the most fundamental question still present a formidable task for analytical methods. Even in the context of classical soliton theories, one of the gems of mathematical physics which culminated with the development of the (inverse) scattering techniques [17][18][19], neither the direct nor the inverse problem generally permit closed-form solutions, and only rare instances are known where the integration can be carried out in an analytic fashion [20][21][22][23]. Indeed, even from a numerical standpoint, no effective framework for computing equilibrium averages of dynamical (or even static) correlation functions is available at this time.
To circumnavigate some of these inherent limitations it is fruitful to attempt a slightly different approach. To better understand certain peculiar features of integrable dynamical systems subject to non-trivial global symmetry constraints, we confine ourselves in this paper to a certain class of classical models in a discrete space-time geometry by following the spirit of a preceding work [24]. Our aim is to explore the possibility of realizing simple exactly solvable symplectic circuits which possess conserved non-Abelian currents. While sacrificing timetranslational symmetry may at first glance seem an unnecessary hindrance, we wish to argue nonetheless that dynamical systems in discrete time offer certain advantages over Hamiltonian models that can be fruitfully employed in various physics applications. An example of this are simple deterministic cellular automata studied recently in [25,26,[26][27][28] which permit one to obtain very explicit results for dynamical correlation functions. In this work, we describe a simple procedure to obtain a class of many-body propagators composed of two-body sympletic maps which governs a discrete space-time evolution of interacting matrix-valued degrees of freedom. This is accomplished in a systematic manner, employing the methods of algebraic geometry and the notion of Lax representation [2] which ensures integrability of the model from the outset. An explicit integration scheme we managed to obtain provides a versatile numerical tool which facilitates efficient numerical simulation of statistical ensembles.
This work has been largely inspired by the recent discovery of superdiffusive magnetization transport in the isotropic Heisenberg spin-1/2 chain [64][65][66], subsequently scrutinized in a number of papers [41,42,67] which collectively piled up a convincing numerical evidence for the Kardar-Parisi-Zhang (KPZ) type universality [68] (see also [69,70] and [71][72][73]). It is remarkable that the same phenomenon is already visible at the classical level, namely in the integrable classical spin chains symmetric under global SO (3) rotations [24,74,75]. In spite of a phenomenological picture based on an effective noisy Burger's equation proposed recently in [41] and further refined in [42], a complete and quantitative understanding of this curious phenomen is still lacking at the moment. More specifically, aside from partial analytical [37,39,76] and numerical evidence [77], it is not very clear what is the precise role of global non-Abelian symmetry and, particularly, if higher-rank symmetries could potentially alter this picture and possibly unveil new types of transport laws. With these questions in mind, we design a class of integrable models whose degrees of freedom are matrix fields which take values on certain compact manifolds. To gain better insight into anomalous nature of spin/charge dynamics in models invariant under the action of non-Abelian Lie groups, we carry out a detailed numerical study of charge transport in maximum entropy states. Our results indicate, quite remarkably, that KPZ scaling is ubiquitous phenomenon independent of the symmetry structure of the local matrix manifold.

Integrable matrix models 2.1 Discrete zero-curvature condition
To set the stage, we shall first introduce the setting. We consider a discrete space-time in the form of a a two-dimensional square lattice. Throughout the paper we adopt the convention that the time flows in the vertical direction and the spatial axis is oriented horizontally towards the right. To each site of the space-time lattice we attach a physical variable. A precise specification of physical degrees of freedom, alongside the associated classical phase space, will be a subject of Section 2.3. By rotating the space-time lattice by 45 • degrees we introduce the light-cone lattice and assign to each of its vertices (nodes) (n, m) ∈ Z 2 an auxiliary variable φ n,m . Physical variables M t are situated on the nodes of the space-time lattice ( , t) ∈ Z 2 (at the midpoints of edges of the light-cone lattice) and we label them as M t=n−m =n+m+1 (resp. M t=n−m+1 =n+m+1 ) when + t is odd (resp. even). We furthermore impose the periodic boundary condition in the space direction, that is M t ≡ M t +L , assuming the system length L to be even.
The outlined construction rests on the notion of a linear transport problem for the auxiliary field, see e.g. references [2,5,18,19,78]. Parallel transport along the light-cone directions (i.e. characteristics x ± t = const) reads φ n+1,m = L (+) n,m (λ)φ n,m , φ n,m+1 = L (−) n,m (µ)φ n,m , where a pair of 'matrix propagators' L (±) , called the Lax pair, represent certain matrix functions of physical variables which additionally depend analytically on the so-called spectral parameters λ and µ. The consistency requirement for the above auxiliary linear problem is that the shifts of the light-cone coordinates commute, meaning that going from φ n,m to φ n+1,m+1 does not depend on the order of the light-cone propagators. This condition can be neatly encapsulated by the discrete zero-curvature property (cf. refs. [79,80] the light-cone lattice Here each of the light-cone Lax matrices L ± depends a local 'edge variable' M , F is a constant invertible 'twisting matrix', and primed variables M are a shorthand notation for M shifted by one unit in the time direction (as depicted in Figure 1). Importantly, the discrete curvature is everywhere zero if and only if the updated variables M 1,2 are appropriately linked to M 1,2 . Another interpretation is to think of the flatness condition as a specification of the dynamical propagator, i.e. a local two-body map (M 1 , M 2 ) → (M 1 , M 2 ) over a spatially adjacent pair of sites (1,2). In Section 2.4.1 we explain how Eq. (2.2) gives rise to integrability of the model. Obtaining and classifying all physically admissible solutions to Eq. (2.2) is likely a difficult task and we shall not undertake it in this work. With a more modest goal in mind, we will attempt to find first the simplest solutions by making the following restrictions: 1. We set both light-cone Lax operators to be equal, L (+) ≡ L (−) .
2. Lax matrix L(λ; M ) is assumed to be a linear function of the spectral parameter λ.
3. Lax matrix L(λ; M ) is assumed to have a linear dependence on the matrix variable M .
We shall interpret a local physical variable M as a classical matrix field which takes values in GL(N ; C) or a submanifold thereof. The third requirement can then be naturally satisfied (without loss of generality) by imposing a non-linear constraint We make the following ansatz for the Lax matrix complying with (1.-3.), and proceed to look for the solutions of the discrete zero-curvature condition of the form It turns out that this matrix equation admits a unique non-trivial solution of the difference type, i.e. that there exist a map (M 1 , M 2 ) → (M 1 , M 2 ) depending solely on the difference of the two spectral parameters µ−λ. As will be demonstrated, the difference condition naturally implies a dynamical conservation law which in the absence of twist (F = 1) implies a global conservation law M t = const when extended to the space-time lattice, see Figure 2. Figure 2: Fabric of discrete space-time: the physical space-time lattice, comprising matrix degrees of freedom M t (blue circles), coexisting with the light-cone square lattice depicted by a tilted checkerboard. A two-body symplectic map Φ τ (red square), which is attached to the middle of each tile, provides the time-propagator for every pair of adjacent physical variables.

Dynamical map
The solution to the zero-curvature condition (2.5), supplemented with nonlinear constraint (2.3), admits a unique solution of the 'difference form' as one-parameter family of symplectic maps, Φ τ : representing diffeomorphisms on the product of two manifolds of involutory matrices. An explicit realization of Φ τ is an adjoint mapping generated by an invertible 1 matrix where the 'twist field' F can be any constant invertible GL(N ; C) matrix. For the proof with a derivation we refer the reader to Appendix A. Note that the map is well defined even for arbitrary complex τ , while we require τ to be real in this paper in order to allow for its interpretation as a Trotter time step of a Hamiltonian flow. The above mapping plays a role of the two-site time-propagator (with time-step τ ) which provides the basic building block of the many-body symplectic circuit (shown in Figure 2) which, moreover, manifestly preserves the non-linear constraint (2.3). The above construction is arguably the simplest integrable many-body dynamical system of non-commuting variables in discrete space-time.
Many-body propagator. The two-body propagator defined in Eq. (2.7) constitutes a basic element of the dynamical map Φ full τ : M L → M L defined on the entire phase space M L = M ×L 1 , denoting the Cartesian product of L copies of M 1 . By virtue of the light-cone structure, the dynamics decomposes into odd and even timesteps, respectively, as depicted in Figure 2. With the usual embedding prescription, where I : M 1 → M 1 designates a local unit function I(M ) ≡ M , the full propagator for a double time step t → t + 2 can be split as, with the odd/even propagators further factorizing as In this view, the full map Φ full τ can be perceived as an 'integrable Trotterrization', closely resembling the integrable Trotterization of the quantum Heisenberg model obtained in [81,82]. Indeed, we will shortly demonstrate in Section 2.6 below that Eqs. (2.8) actually correspond to complete space-time discretizations of (nonrelativistic) σ-models with compact Lie group 1 Non-degeneracy of Sτ for |τ | > 0 follows from showing Det(Sτ )Det(S † τ ) > 0, which is a consequence of hermiticity of M1,2, commutativity [M1M2, M2M1] = 0, M1,2 2 = 1, and sub-multiplicativity of the operator norm.
cosets as their local target spaces. These include, in particular, the higher-rank analogues of the Landau-Lifshitz magnets [18,83] pertaining to CP n manifolds, cf. [84]. This class of classical field theories indeed naturally emerges as the semi-classical limit of integrable quantum chains of locally interacting 'spins' which exhibit manifest symmetry under SU (N ), introduce a long while ago in [85,86]. In Appendix E we show that the effective classical action governing the long-wavelength modes above a ferromagnetic vacuum yields precisely the continuous counterparts of our matrix models.

Phase space and invariant measures
In this section we proceed by identifying the admissible phase spaces for the dynamical map and describe their formal properties. for a general introduction to differential geometry and symplectic manifolds we refer the reader to one of the standard texts, e.g. [87,88].
In classifying the phase space, the non-linear constraint (2.3) plays a pivotal role. Moreover, it will be key to us to restrict ourselves to compact smooth manifolds where the notion of a normalizable invariant measure is well defined. This leaves us with the following list of compact Lie groups: (i) unitary groups U (N ), (ii) orthogonal groups O(N ), and (iii) compact symplectic groups USp(2N ). 2 Complex Grassmannians. We begin be examining the unitary Lie groups G = U (N ), assuming N ≥ 2. The first thing to notice is that, by virtue of the involutory property (2.3), the dynamical matrix variables do not bijectively correspond to the elements of G but instead lie on a submanifold spanned by N -dimensional hermitian matrices prescribed by (2.14) This defines the so-called complex Grassmannian manifold, the set of k-dimensional complex planes embedded in C N which pass through the origin, corresponding to the eigenspaces of M with eigenvalue −1. Specifically, by prescribing a diagonal signature matrix, N ) can be naturally identified with the adjoint orbits of Σ (k,N ) under the action of group G, that is any M ∈ M can be obtained as (2.16) Here we emphasize that a matrix M ∈ M 1 , as given by Eq. (2.16), does not correspond to a unique group element g, the reason being that Σ (k,N ) is invariant under conjugation with unitary matrices of the form In this view, Gr C (k, N ) are homogeneous spaces, i.e. cosets of the isometry group G by the stability group H, , (2.17) where ∼ = stands for diffeomorphic equivalence. 3 The group G acts transitively on each component of Gr C (k, N ) by virtue of Eq. (2.16), with the subgroup H being the G-stabilizer of Σ (k,N ) . In other words, the group manifold foliates into equivalence classes under the action of H, each of which is then understood as elements of the coset space G/H. Grassmannian manifolds include, as a special case, complex projective spaces Gr C (1, N ) ∼ = CP N −1 , namely complex lines passing through the origin of a complex Euclidean space. Due to equivalence Gr C (k, N ) ∼ = Gr C (N − k, N ) we shall subsequently assume, without loss of generality, that k ≤ N/2 . Matrices which satisfy the involutory property (2.3) can be alternatively realized in terms of rank-k projectors P projecting onto the eigenspace with eigenvalue −1, A brief inspection shows that real Grassmannians Gr R (k, N ) are not preserved under the dynamical map (2.8) and hence will not be considered further in this work.
Lagrangian Grassmannians. We finally consider complex matrices which preserve the symplectic unit, namely the symplectic group Sp(2N ; C). Although the latter is not compact, its restriction to the unitary subgroup, that is the intersection of Sp(2N ; C) with SU (2N ), is a simply-connected compact Lie group USp(2N ) of 2N -dimensional complex matrices (2.19) called the unitary symplectic group. Here J denotes the standard symplectic unit The adjoint USp(2N ) orbits of the antisymplectic signature 4 In Appendix A we demonstrate that the dynamics (2.8) preserves the Lagrangian submanifold with the anti-symplectic signature, hence Lagrangian Grasmannians again also constitute an admissible phase space M 1 ∼ = L(N ).

Affine parametrization
To specify a Grassmannian manifold Gr C (k, N ) one has to supply k linearly independent complex vectors of dimension N . Let us suppose these are stored as columns of a complex N ×k matrix Ψ. Matrix elements of Ψ are referred to as homogeneous coordinates of Gr C (k, N ). It is crucial to recognize here that the choice of basis vectors is not unique as one enjoys the freedom of performing linear transformations Ψ → ΨA with any invertible k-dimensional matrix A. Given that Grassmannians are identified with equivalence classes g H, a description in terms of homogeneous coordinates involves (in general non-Abelian) gauge freedom. By exploiting this gauge redundancy, we can always pick and choose an element in each equivalence class to bring the coordinate matrix into a 'canonical form' Grassmannian manifolds Gr C (k, N ) therefore corresponds to complex manifolds of real dimension dim Gr C (k, N ) = dim U (N ) − dim (U (k) × U (N − k)) = 2k(N − k) parametrized locally by affine coordinates z i,a . By parametrizing the group element in the form [89][90][91] a rank-k projector P (Z) assumes the block structure (2.27) One shortcoming of such an explicit parametrization is that it does not provide a global parametrization of M 1 . Indeed, one needs in total N k coordinate charts to cover the entire phase space, obtained by all possible distributions of −1 in the diagonal signature matrix Σ.

Symplectic structure
Complex Grassmannians Gr C (k, N ) are symplectic manifolds. They are indeed Kähler manifolds, meaning that they possess compatible Riemannian and symplectic structures. In this section we give a succinct review of the basic notions which we shall subsequently use throughout the rest of the paper.

SciPost Physics
Submission Sympletic form. The local phase space M (k,N ) 1 ∼ = Gr C (k, N ) of a matrix model is a smooth manifold endowed with a symplectic 2-from ω K which is closed, dω K = 0, and non-degenerate, Det(ω K ) = 0. Expressed in terms of local affine coordinates z i,a and their complex conjugates z i,a (with i = 1, 2, . . . , N − k and a = 1, 2, . . . k), the Kähler form ω K reads compactly The Kähler form can be expressed in terms of the Riemann metric tensor, 29) and, employing the vectorized matrix coordinate can be written as ω K = 1 2i dZ † ∧ η(Z)dZ. In the special case of complex projective spaces CP n this reduces to the well-known Fubini-Study metric η FS , which can be obtained from the Kähler potential An alternative way of introducing the symplectic structure is to exploit the algebraic structure. This is not only advantageous from the practical standpoint, but also does not require any particular coordinatization of M (2.32) We can again avoid making any reference to an explicit coordinate system and use the fact that there is a natural action of the group G = SU (N ) on M 1 given by conjugation M → g M g † , where g ∈ G is a group element of the form g = exp (−i a θ a X a ), with θ a ∈ R and X a ∈ g being traceless hermitian matrices generating the Lie algebra g = su(N ). Fixing the basis {X a }, a = 1, 2, . . . , dim g = N 2 − 1, with normalization where abc is the appropriate tensor of structure constants. Now the matrix-valued vector fields V X (M ) can be viewed as an infinitesimal action of group G at M ∈ M 1 , that is Moment maps. The Lie algebra structure on the space of vector fields realized by the commutator induces a Lie algebra structure on the space of functions provided by the Poisson bracket on the phase space M 1 . A mapping from a Lie algebra to functions on classical phase spaces is given by the moment map, formally obtained by contracting the symplectic form with the vector field df X = ι V X ω. (2.36) Contracting the symplectic form using dM (V X ) = V X (M ), In any local coordinate chart, the Poisson bracket can be expressed through the inverse of the Riemann metric To avoid an explicit coordinate description we instead utilize the moment maps associated to the action of g and accordingly define the Poisson bracket through the contraction of ω, This readily implies the following relation for the moment maps, yielding the Lie-Poisson algebra, Using furthermore that we deduce the su(N ) Lie-Poisson algebra for the hermitian components of M , At this point we wish to stress once again that, by virtue of M 2 = 1, not all matrix elements (components) M a appear as independent fields. In fact, a nonlinear constraint amounts to fixing a symplectic leaf (i.e. Casimir invariants), rendering the symplectic form ω non-degenerate. The Lie-Poisson bracket (2.44) on a local phase space M 1 can be more conveniently reformulated in an equivalent matrix form where Π is the permutation (transposition, or swap) matrix over C Finally, promoting the Lie-Poisson algebra to Lax matrices, the linear bracket (2.46) can be presented as the Sklyanin's fundamental quadratic bracket where the intertwiner r(λ, λ ) is the so-called classical r-matrix [12,18,19].
Hamiltonian field. The Hamiltonian action associated with a vector field V H (H ∈ g) induces the following dynamics of the moment maps At the level of matrix variables M , one can accordingly deduce the 'Heisenberg equation of motion' d dt and using the projector representation of the moment map, the Hamiltonian equations of motion can also be given in the affine coordinates (cf. Eq. (E. 7) in Section E.1) along with the complex-conjugate counterpart. These can be recast compactly in the form of a matrix Ricatti equation [89] d dt

Separable invariant measure
Since M (k,N ) 1 are homogeneous spaces, they admit a G-invariant measure inherited from the invariant Haar measure of the unitary group G. This measure is nothing but the normalized invariant symplectic volume, defined via the highest exterior product of the Kähler form ω K with itself In terms of the Riemann metric tensor we therefore have The determinant of the metric tensor can be expressed in terms of the affine matrix coordinate (2.57) Liouville measure. The Liouville measure specified by density ρ (k,N ) refers to an unbiased (flat, or uniform) G-invariant probability measure on M (k,N ) 1 given by the normalized symplectic volume, The Liouville volume N (k, N ) can be computed in an indirect manner by exploiting the coset structure of Gr C (k, N ). The volume of the unitary group U (n) can easily found via the isomorphism U (n)/U (n − 1) ∼ = S 2n−1 , where the volumes of hypersphere are known to be Vol(S 2n−1 ) = 2π n /(n − 1)!, whence Using the coset structure we thus have [92] In physics applications, the above measure can be naturally related to the maximumentropy (or infinite-temperature) equilibrium ensemble of a local matrix degree of freedom, namely a distinguished measure which maximizes the Shannon/Gibbs entropy (2.61) The Liouville measure over the many-body (product) phase space M (k,N ) L is then given by a product (separable) flat measure To establish that the Liouville measure (2.62) is invariant under the time evolution generated by Φ full τ , it suffices to verify that the symplectic form ω (and hence the volume element dΩ) is preserved under the action of two-body propagator Φ τ given by Eqs. (2.8). This amounts to demonstrating that Φ τ preserves the Poisson bracket, as shown explicitly in Appendix B.
Grand-canonical measure. The Liouville measure on M (k,N ) 1 admits a multi-parameter extension with the normalization factor interpreted as the grand-canonical partition function. This is formally a push-forward of the Liouville measure by the moment map known in the mathematical literature as an equivariant measure. The class of grand-canonical measures (2.63) solves the constrained variational problem of entropy maximization (cf. Eq. (2.61)) with prescribed Lagrange multipliers µ a ∈ R.
Again we can build a product measure over M L from the grand-canonical measures over single sites. Since any two adjacent local phase spaces M 1 × M 1 in the Cartesian product M L are acted on by the group G in a Hamiltonian fashion and the action of G is diagonal, an equivariant measure ρ (k,N ) L on M L is, provided F = 1, also preserved under the action of time-propagator Φ τ as an immediate consequence of the fundamental conservation law (2.6).
We can assume, with no loss of generality, that the grand-canonical measure is characterized only by the maximal torus of G, namely that the exponent in Eq. (2.63) is an element X a ∈ g 0 of the maximal Abelian (Cartan) subalgebra g 0 of g. The phase-space averages of the Cartan charge densities are given by In performing the phase-space integration over an Abelian equivariant measure ρ(M ; {µ a }) explicit integration can be avoided thanks to the localization theorem due to Duistermaat and Heckman [93] (see also [94,95]). The latter is essentially a statement about the exactness of the saddle-point approximation: an integral of the exponent of the torus action on a compact phase space localizes at its critical points. Using the fact that the torus action on M 1 is governed by a linear matrix equation, where A and D are blocks of Hamiltonian matrix (2.51), with H = a µ a X a . Assuming for definiteness that the critical points are all isolated (i.e. non-degeneracy A a,a = D i,i for all a = 1, . . . , k and i = 1, . . . , N − k), the stationary points (d/dt)Z = 0 in every coordinate chart are located precisely at the origin Z = 0, with the Hessian matrix Denoting χ α = H αα and applying the Duistermaat-Heckman formula one finds (see e.g. [90]) where the summation over ordered sets σ = {σ 1 < σ 2 < . . . , σ k } covers all N k coordinate patches (i.e. all possible redistributions of −1's in the signature Σ (k,N ) ) andā runs over the complementary set of indicesσ = {1, 2, . . . , N } \ σ.

Isospectrality
In this section we discuss certain integrability aspects of Eqs. (2.8) which arise as direct consequences of the zero-curvature representation (cf. Section 2.1) on the discrete two-dimensional light-cone lattice. The flatness property of the Lax connection implies that parallel transport of the auxiliary variable from one point on the light-cone lattice to another does not depend on the path (Wilson line) between the two. This means that all contractible closed paths on the light-cone lattice are trivial. On the other hand, a discrete holonomy corresponding to a non-contractible path which wraps once around the system (with periodic boundary conditions) is non-trivial and provides an analytic family of matrix-valued functions on the phase space M L called the (staggered) monodromy matrix where we have adopted the right-to-left path-ordering convention. Owing to the zero-curvature condition, monodromies evaluated at a pair of time arguments of equal parity are related to one another by a similarity transformation, implying that their eigenvalues or any spectral invariants are conserved under the time evolution. Since M τ (λ|{M }) depends analytically on the spectral parameter λ, it provides a generating functional for an extensive number (i.e. O(L)) of constants of motion (conservation laws), namely functionally independent phasespace functions preserved under the evolution. This hallmark property of Lax integrability is referred to as isospectrality [19,79]. Establishing integrability in the Liouville-Arnol'd sense amounts to showing in addition that these conservation laws are also mutually in involution.
The trace of the monodromy matrix defines the transfer map T τ (λ) : M L → C, constituting a family of analytic phase-space functions on M L mutually in involution 9 where τ ∈ R is a fixed parameter. By sequential application of the local zero curvature relation (2.5) on, subsequently, even and odd pairs of Lax operators along a fixed horizontal sawtooth on the light-cone lattice, one obtains time conservation of the transfer map T τ (λ) then generates, via logarithmic differentiation, a family of local conservation laws, at least for signatures with rank k = 1. This is shown explicitly in Appendix D.

Space-time self-duality
Dual propagator. The two-body propagator defined in Eq. (2.7) realizes the time propagator of the matrix model, namely it propagates a pair of adjacent matrix variables by one unit step in the time direction. The discrete zero-curvature property (2.5) however also permits to define the dual propagator Φ d τ : M 1 × M 1 → M 1 × M 1 , the spatial analogue of a two-body map where variables (M 1 , M 1 ) appear as an 'incoming' state and (M 2 , M 2 ) as an 'outgoing' state. See Ref. [96] for a discussion of related concepts in continuum integrable models.
To construct the dual propagator, we consider the linear transport problem for the auxiliary fields in the time direction, noticing that reversing the direction of propagation amounts to inverting the Lax operator (cf. Eq. (2.4)) Starting from Eq. (2.5), acting from the right by the inverse of L(µ; M 1 ) and from the left by the inverse of L(µ, M 2 ), and finally multiplying by twists F −1/2 from both sides, we obtain the dual zero-curvature relation (as depicted Figure 3) A neat trick to solve this equation is to apply a local gauge transformation, which brings us back to the original form of the (twisted) discrete zero-curvature relation (2.5) along with the dual dynamical symmetry law (2.6) The upshot here is that the local two-body spatial propagator in the 'tilde' (gauged) variables exactly coincides with the temporal one: An explicit prescription for the dual propagator Φ d τ , can be found by undoing the gauge transformation (2.75), In conclusion, the dual (i.e. spatial) propagator Φ d τ and the temporal propagator Φ τ (2.7) are, apart from a local gauge transformation of the plaquette, identical maps.

Self-duality.
In the absence of twist (F = 1), the above gauge transformation can be consistently extended to the entire space-time lattice implying that in tilde variables the full spatial dynamics can be expressed in terms of the temporal propagator We shall refer to this property as space-time self-duality.
Two remarks are in order at this point. First, it is worth pointing out that the self-duality property, despite its manifest presence in the fully discrete setting, disappears at the level of the Hamiltonian dynamics which emerges in the continuous time limit (cf. Eq. (2.95)). This is caused by the fact that space and time coordinates no longer appear on equal footing in the continuous time limit. Indeed, in deriving the continuum limit one only retains smooth variations of the classical field configurations, which is clearly in conflict with the staggered nature of the local gauge transformation (2.75).
The self-duality property also somewhat insinuates that dynamics in the matrix model may be subjected to certain stringent restrictions. It may appear, at the first glance at least, that the temporal and spatial dynamics in the model should be effectively interchangeable. This is not at all what happens. For instance, an uncorrelated time-invariant initial state after being locally quenched undergoes a non-trivial time-evolution resulting in a strongly correlated 'time state' [28]. Furthermore, correlations do not only propagate with the unit speed as e.g. in the dual-unitary models [97]. In Section 3 we carry out numerical simulations to explicitly demonstrate this fact. To recover the canonical form we additionally apply a gauge transformation on the matrix variables (bottom panel), thus establishing a duality between the temporal and spatial propagators.

Yang-Baxter relation
Another important manifestation of integrability (cf. the zero-curvature property (2.5)) is that the two-body propagator is also a Yang-Baxter map, R λ : whereΦ λ denotes the untwisted elementary propagator and Π is the permutation map on M 1 × M 1 . By embedding the maps into a triple Cartesian product M 1 × M 1 × M 1 , we find that R λ satisfies the set-theoretic Yang-Baxter relation whereas the untwisted propagator accordingly satisfies the associated braid relation Let us briefly elucidate the origin of the Yang-Baxter map (see [79,[98][99][100], or [101,102] for more recent accounts which discuss its connection to quasitriangular Hopf algebras [103,104]). To this end it is convenient to think of the discrete zero-curvature condition (2.5) as a refactorization problem [105] for a pair of Lax matrices where ζ j ∈ C are arbitrary shift parameters. The Yang-Baxter map R λ provides a map- which is a unique solution to the matrix re-factorization problem. The set-theoretic (functional) Yang-Baxter property is a statement about equivalence of two different intertwining protocols; applying the left-and right-hand sides of Eq.

Symplectic generator
Having shown that Eq. (2.7) provides a symplectic transformation, we can alternatively realize it in the Hamiltonian form. To this end we define is not simply an integrable lattice Hamiltonian obtained in the τ → 0 limit (derived below in Section 2.6). It is also important to stress that in the absence of time-translational symmetry the notion of energy is not meaningful. In this respect, the generating function H is merely an auxiliary concept which need not necessarily be a real function in general. The only physical requirement, besides generating the symplectic map (2.8), is that in the limit of continuous time the generator H indeed yields a real (integrable) Hamiltonian function.
For definiteness we confine ourselves here to the untwisted case and set F = 1. The generator H is in general a functional of trace invariants of S 0 = M 1 + M 2 , that is scalar functions s m = Tr(S m 0 ). Taking into account that s 0 = N , s 1 = 2(N − 2k), and that all odd invariants s 2m+1 are proportional to s 1 as a consequence of cyclicity of the trace and M 2 1,2 = 1, is thus sufficient to only retain only s 2m for m ∈ N. We have succeeded in deriving a system of PDEs that determines the generator H wheres is defined through p 1,3 (s) = 0, which contains an extra purely imaginary part. The latter nonetheless vanishes in the limit τ → 0, which we expect to be a general feature of symplectic generators in odd-dimensional cases.
With aid of the Vieta's formulas and the Jacobi identity we moreover infer the local Hamiltonian in the continuous time limit, again for even N :

Semi-discrete and continuum limits
To better elucidate the physical meaning of our matrix models it it is instructive to also inspect their time-continuous and field-theoretical limits. We consider first the limit τ → 0, where the symplectic map Φ τ 'smoothens out' into an integrable Hamiltonian flow governed by a lattice Hamiltonian H lattice . For this purpose we parametrize the twist field as F = exp(−iτ B/2), with B ∈ g, and expand Eq. (2.7) to the lowest order in τ , yielding the differential-difference equation The Poisson bracket for this field theory is found by taking the continuum limit of the linear Poisson bracket (2.46), Lax equations. The auxiliary linear transport problem for the auxiliary field φ (t) in discrete space and continuous time takes the form where the spatial propagator L is the Lax matrix (2.4) inherited from the light-cone lattice and V denotes the temporal component which we determine below. The compatibility condition for Eqs. (2.98) takes the form of a semi-discrete zero-curvature condition, (2.100) To obtain the Lax connection for the continuum counterpart, we reintroduce the lattice spacing parameter ∆ and expand Eq. (2.99) to the second order O(∆ 2 ). In this way, we obtain the auxiliary linear transport problem associated to a differentiable manifold, satisfying the zero-curvature compatibility condition with connection components Example. For a brief illustration, we consider the simplest example of a 2-sphere M 1 = S 2 ∼ = CP 1 . In this case, it is useful to represent the matrix field variable M in terms of a unit vector (spin) field S ∈ S 2 (S · S = 1) in R 3 , M = S · σ, where σ = (σ x , σ y , σ z ) T is a vector of Pauli matrices. The symplectic map Φ τ for this case has been studied in the previous work [24] Φ τ (S 1 , To reproduce the semi-discrete and continuum limit of Eq. (2.104) we expand the inverses in Eq. (2.93), yielding where we put B = B · σ/2. This is nothing but the integrable lattice discretization of the SO(3)-symmetric Heisenberg ferromagnet (isotropic lattice Landau-Lifshitz model [106,107]) in a homogeneous field B. The equation of motion is generated by the logarithmic interaction of the form [18]

Charge transport and KPZ superuniversality
The remainder of the paper is devoted to numerical study of transport properties of integrable matrix models in canonical equilibrium ensembles. Since our aim is mainly to address the central questions outlined in the introduction, we shall focus entirely on the Noether charges where certain anomalous features are to be expected. In particular, we refer to previous studies of magnetization transport in the isotropic Landau-Lifshitz (Heisenberg) magnet (with the spin-field belonging to the coset G/H = S 2 ) which uncovered superdiffusive transport of the KPZ universality class both in the quantum and classical setting [24,[64][65][66][67]74,75]. The aim of the subsequent analysis is to systematically analyze the role of isometry and isotropy groups G = SU (N ) and H = S(U (k) × U (N − k)), respectively. We shall also consider a distinct case of symplectic symmetry with G = USp(2N ) and H = U (N ). The Noether charge is a G-valued dynamical observable whose local densities are provided by the moment map f X (M t ) = Tr(X M t ). (3.1) For the basis of the Noether charges we will subsequently use notation q a (t) ≡ f X a (M t ).
Exact computation of time-dependent correlation functions in equilibrium states lies beyond the capabilities of available analytic techniques. We thus have to rely fully on numerical simulations. The main object of study in our simulations are connected spatio-temporal autocorrelation functions of charge densities q a (t), where the 'equilibrium expectation value' · pertains to averaging with respect to a uniform Liouville measure on M L , an analogues of the canonical Gibbs state at 'infinite temperature', invariant under unit time and space shifts t → t + 1 and → + 1, respectively. Indeed, since G acts transitively on Gr C (k, N ), the G-invariant measure on Grassmannian manifolds is naturally inherited from the invariant (Haar) measure on G. In this view, we can first sample uniformly over the group G (see e.g. [108]) and then generate the invariant distribution on Gr C (k, N ) through the mapping M = g Σ (k,N ) g † , see Eq. (2.16).
We have numerically computed the dynamical correlator defined in Eq. (3.2) using the following scheme. First, we generated N s initial matrix ensembles E ≡ {M t=0 } Ns α=1 by drawing each sample set from the Liouville probability measure ρ (k,N ) . Next, we computed the connected longitudinal dynamical correlators with the following prescription which can be efficiently performed using the convolution theorem. The maximal time of simulation t max can be adjusted so as to eliminate any spurious effect due to periodic boundary conditions. 12 To smear out the even-odd effect of staggering (see Figure 2), it is better to compute the autocorrelation function of the Noether charges by averaging over adjacent pairs of variables, q := 1 2 (q + q +1 ), corresponding the 'smoothened' correlation function Lastly, by virtue of the global G-invariance we are allowed to average over all the components a = 1, 2, . . . , dim g.

Uniform equilibrium states
In Figure 5 we display the time-dependent correlation functions (3.2) in the absence of twist field (F = 1) for the few smallest dimensions N ∈ {3, 4, 5} and all inequivalent signature Assuming an algebraic decay at large times, we first extract the dynamical exponent z from the numerical data. We find, uniformly for all the instances with k ≤ N/2 and N = 2, 3, 4, . . ., excellent agreement with the Kardar-Parisi-Zhang superdiffusive universal algebraic exponent z KPZ = 3/2, cf. Figure 6.
To further corroborate the presence of the KPZ physics, we proceed with the extraction of the scaled dynamical structure factor C q (ξ, t) = t 1/z C q ( , t), ξ := t −1/z . (3.6) In Figure 7 we display the stationary cross sections which are expected to collapse onto a universal scaling function g P S tabulated in [109], where λ B ∈ R is the Burger's field coupling constant and A is the amplitude. Agreement with the universal KPZ profile is very solid and there are clearly visible systematic deviations from the Gaussian form which is characteristic of normal diffusion, see Figure 7. The non-Gaussian behaviour of the scaling function g PS is most pronounced in the tails, i.e. at large values of the scaling variable ξ. The extracted numerical values of constants λ B and A are reported in Table 1.
For completeness we include the numerical analysis of charge transport in an integrable matrix model on a Lagrangian Grassmannian. Since these are sub-manifolds of complex Grassmannians, they have to be considered independently. We shall only consider here the simplest instance L(2) ∼ = USp(2; C)/U (2). Numerical data shown in Figure 8

Magnetic field
To incorporate an external SU (N )-magnetic field we set the twisting element to where B is a fixed Hermitian matrix of the form B = h a n a X a , with field strength h and (unit) polarization vector n. The addition of a field causes the following dynamical effect: the dynamical correlations pertaining to the distinguished Noether charge aligned with the polarization direction is unaffected by the field, whereas the correlation function for all the remaining charges exhibit a super-diffusive KPZ spreading modulated by a periodic precessional motion, see Figure 9.

Inhomogeneous phase space
We finally mention an interesting possibility of introducing an integrable matrix model on an inhomogeneous phase space which admits the general structure which is still compatible with integrability of the many-body dynamics. This is a corollary of the fact that Φ τ acts as a conjugation in G × G which preserves the the total signature by swapping the signature of two adjacent incidence matrices M and M , and there is no issue with 'scattering' degrees of freedom from different adjoint orbits. As a consequence, the total signature Σ (k ,N ) is conserved under time evolution.

Staggered phase space.
As an illustration of the above construction we consider a special case of the staggered phase space with an alternating sequence of inequivalent phase spaces of rank k = 1 and k = 2, specializing to the lowest-dimensional instance N = 4. As shown in Figure 10  . Amplitudes A i were read off from the horizontal axis intercepts of the equal-space correlator C(0, t) shown in Figure 6. The last line pertains to the unitary symplectic case (see Figure 8).
appear to converge towards the KPZ scaling function. In fact, we find an asymmetric profile with discernible deviations in the left tail which do not seem to come from finite-time effects. Given that the model with a staggered phase space has an intrinsic chiral structure, we have also tried a two-sided fit, that is fitting the KPZ scaling function for the left and right movers separately. Doing this however did not appreciably improve upon the fit in Figure 10. We postpone a more detail analysis of this exceptional scenario for future work.

Discussion and conclusion
We have introduced a novel family of integrable models of matrix-valued classical fields propagating on a discrete space-time lattice, and obtained an explicit dynamical system in the form of a classical Floquet circuit composed of elementary two-body symplectic maps. The class of models is distinguished by the presence of a conserved G-invariant Noether currents and (in general non-Abelian) local gauge invariance under a subgroup H. In the absence of external fields, both time and space dynamics can be realized in uniform way, highlighting a particular type of space-time self-duality. Integrability of our models manifest itself through the discrete zero-curvature condition on the light-cone lattice, implying infinitely many conserved quantities in involution and set-theoretic Yang-Baxter relation for the elementary two-body symplectic propagator.
Integrable difference equations have been extensively studied in the mathematical physics literature (see [79,80]), particularly in the context of equations on quadrilateral graphs [99,110,111] which have been classified in the work of Adler, Bobenko and Suris [112]. These includes, among other, the Faddeev-Volkov discretization of the sine-Gordon model [113]. We are nevertheless not aware of any known classification scheme which would comprise the class of models discussed in this work, despite from the viewpoint of conceptual simplicity they can hardly be rivalled. Another known, albeit somewhat less explored, approach to produce integrable difference equations, developed in [114][115][116] and [117][118][119], is via discretization of Hirota derivatives [120,121]. Although with this approach the authors of Ref. [117] succeeded to obtain a particular lattice discretization of the N = 2 isotropic Landau-Lifshitz model, its implicit form makes it less appealing for concrete applications. In more recent works [101,102] several formal connections between quantum Yang-Baxter maps (representing the adjoint action of the universal R-matrix of a quantum group) and their classical limits (and discretetime dynamics) have been outlined, indicating that classical Yang-Baxter maps in a way naturally descend from the associated quantized algebraic structure. In this respect, our results indicated that at the set-theoretic level, the emergent classical Yang-Baxter map does not even depend explicitly on the underlying Lie algebra.
We devised an explicit integration scheme and employed it as an efficient numerical tool to investigate the dynamics of Noether charges in unbiased invariant maximum-entropy states. In analogy with the isotropic Landau-Lifshitz model, we now found robust evidence of the superdiffusive dynamics of the KPZ type irrespectively of the structure of the local phase space, i.e. isometry and isotropy groups of their coset target spaces. In fact, universality of KPZ type extends even to Lagrangian Grassmannians that are linked with unitary symplectic groups. Despite the outline construction does not accommodate for matrix models associated with real orthogonal Lie groups and associated real Grassmannians or exceptional compact groups, these could presumably be included with suitable adaptations.
We nonetheless believe that the following conjecture can be stated: all discrete space-time models built as Floquet circuits from two-body symplectic Yang-Baxter maps (and continuum limits thereof), with dynamical variables living on compact non-abelian symmetric spaces exhibit superdiffusion of the KPZ type in equilibrium states with unbroken symmetry. If that is indeed the case, the observed phenomenon of KPZ physics can be dubbed as superuniversal. This conjecture could be even slightly expanded by adjoining also matrix models on supersymmetric coset spaces invariant under Lie superalgebras (it is known from [37] that the associated integrable quantum chains have singular diffusion constant as well). At the moment we postpone further examination and other related questions to future studies, including a more comprehensive study of charge transport by extending the analysis to grand-canonical states.
Since our models provide integrable Trotterizations of (non-relativistic) coset σ-models with complex Grassmannian manifolds as their target spaces, which emerge as semi-classical limits of integrable symmetric quantum spin chains invariant under global SU (N ) rotation. Since charge transport in these integrable quantum chains is also known to be of anomalous type [37], it appears plausible that our findings translate directly to the quantum setting too. This would suggest that there is a general principle behind an exact quantum-classical correspondence of charge transport, discussed in the scope of the Heisenberg model in Refs. [23,122].
There are several distinct features of our models that were left unexplored but definitely merit further study. On the formal side, developing a fully-fledged inverse scattering formalism by integrating the auxiliary linear problem in discrete space-time would provide a platform to tackle various problems of nonequilibrium statistical mechanics in an analytical fashion. Another pending question is the fate of local conserved quantities beyond the projective models (k = 1). Since for Gr C (k ≥ 2, N ) the monodromy matrices evaluated at the projection points no longer decompose into a sequence of rank-1 projectors, it is not obvious which mechanism (if any) would ensure locality of conserved quantities. Curiously however, in the limit of continuous time the symplectic generator of the elementary propagator yields a strictly local Hamiltonian density governing the time evolution of their integrable lattice counterparts even for generic (i.e. non-projective) models (k ≥ 2). It is not inconceivable that these lattice models involve quasi-local conservation laws, bearing some resemblance to higher-spin commuting transfer matrices in quantum Heisenberg model where the 'shift point' property also ceases to exist [123,124]. Another interesting question which remains open is whether it is possible to take different continuum limits to systematically recover the entire hierarchy of higher Hamiltonian flows in the field-theory limit. While the space-time self-duality property of the symplectic map permits us to study dynamics in the space direction, that is evolutions of time-states (see e.g. [28]), it is less clear at this moment if this has any implications on the structure of spatio-temporal correlations in these matrix models or possibly even on the observed anomalous transport of Noether charges. It has to be stresses that this type of self-duality differs fundamentally from the so-called dualunitarity found recently in the context of quantum circuits [97,125] where correlations are locked to the light-rays. In our models, correlations function with respect to a flat invariant measure fill the entire causal cone in a non-trivial fashion.
From the viewpoint of other physical applications it would be valuable to obtain various integrable deformations such as adding an uniaxial interaction anisotropy. Based on numerical evidence from the anisotropic Landau-Lifshitz model [74], spin transport depends quite intricately on the value of anisotropy, ultimately responsible for the elusive behavior of the gapless phase in the anisoropic Heisenberg XXZ spin-1/2 chain [126,127] (see [124] for a review). To conclude, we wish to mention that the observed superuniversal nature of the KPZ phenomenon in integrable models with non-Abelian Noether currents gives a further (albeit implicit) hint that the notion of hydrodynamic soft modes, employed recently in a phenomenological description of the KPZ phenomenon [41,42], extends beyond the simplest SO(3)-invariant Landau-Lifshitz theory. We postpone the study of these aspects for the future.

Appendices A Two-body propagator
Here we show that the discrete zero curvature condition Equation (A.1) represents equality of two quadratic polynomials in the spectral parameter λ ∈ C. We must thus equate each power in λ. Expanding out the zero-curvature condition (A.1) and dropping the leading λ 2 terms, we find Matching the linear terms in λ immediately implies a two-site conservation law (2.6) while the λ 0 term gives Using Eq. (A.3), the later can be cast as from where it follows Plugging the above back into Eq. (A.4) to eliminate the product of primed variables and writing S τ = M 1 + M 2 + iτ , we arrive at the following explicit form In this way, we have established uniqueness of time propagator Φ τ . The right-hand side can be brought into the final form (2.8) by a straightforward computation. Gr C (k, N ). It can now be readily demonstrated that the mapping (2.8) preserves Grassmannian submanifolds Gr C (k, N ) = SU (N )/S(U (k) × U (N − k)) of the unitary group SU (N ) if τ ∈ R, without explicitly specifying N and k. Considering two hermitian matrices M 1,2 and a unitary twisting matrix F , the hermitianconjugate counterpart of

Closure on
To establish M 1 = (M 1 ) † , is remains to verify explicitly that S −τ S τ M 2 = M 2 S −τ S τ , which is a matter of a short calculation. Preservation of the involutory property under Φ τ is manifest from Eq. (A.9).
Closure on L(N ). Let M 1,2 now be two unitary antisymplectic complex matrices of dimension 2N with the involutory property, and F a symplectic or antisymplectic matrix, namely

B Symplectic properties
Here we present a direct proof that Eq. (2.8) provides a symplectic map (i.e. symplectomorphism) on the product of two Grassmannians manifolds M 1 × M 1 and thus also conserves the product Liouville measure. To establish this, it is enough to demonstrate that the Poisson bracket is preserved under the time-evolution. It proves convenient to carry out this computation at the level of the Sklyanin bracket.
First we treat the bracket operating on a pair of variables on the same space M 1 . Separating out the twist dependence, and subsequently expanding everything out using the Leibniz derivation rule, a tedious calculation yields The obtained expression be further simplified with the aid of general matrix identities valid for any set of dummy N -dimensional matrices {A, B, C, D}. 14 Using these, Eq. (B.2) can be brought into the form in agreement with Eq. (2.45). To establish ultra-locality of the Poisson bracket (cf. Eq. (2.46)) we need to additionally verify that This can once again be confirmed with an explicit but lengthy computation, which, after simplification, indeed evaluates to zero.

B.1 Symplectic generator
The aim of this section is to reformulate the two-body symplectic map ( meaning that only the even invariants s 2k can be non-trivial functions of the dynamical variables.
By application of Leibniz's rule, the equations of motion are therefore put in the form where momentarily no upper limit in the summation has been imposed. The vital piece of the derivation is thus to compute the Poisson brackets {M 1,2 , s 2k }, which can be achieved by using partial traces. Exploiting a useful identity we proceed by evaluating the Poisson bracket This allows us to deduce the equations of motion This expression manifestly demonstrates that the sum M 1 + M 2 is in fact conserved in time, which permits us to write the Heisenberg equation of motion with the solution  Evaluating now the solution at t = τ we deduce a neat time-translation property which rewards us with a remarkably simple expression for F By expanding the logarithm into a series, multiplying by S m 0 and taking the trace, we find an infinite system of linear partial differential equations The solution to this system provides us with the symplectic generator H , uniquely up to additive constants. Indeed, the solution is guaranteed to exist by symplecticity of the timepropagator Φ τ , and a brief inspection shows that equations at different values of n are also independent.
The infinite system (B.22)) can be further reduced to a finite closed system of differential equations by performing a resummation of the matrix invariants s k . This can achieved by means of the Cayley-Hamilton theorem which states that every N -dimensional matrix A satisfies its own characteristic polynomial, Coefficients c j are provided by traces of powers of A, a j ≡ Tr A j , with B j denoting the exponential Bell polynomials. The key observation here is that a matrix of dimension N possesses only N independent scalar invariants (for instance a j up to j = N ). Specifically, we define the Cayley-Hamilton polynomial corresponding to the signature Σ (k,N ) as We can accordingly proceed by solving the following truncated system of partial differential equations We shall not attempt to find its general solution here, but instead rather consider the few simplest instances for small matrix dimensions N . This is in agreement with the expression found previously in [24]. In the τ → 0 limit, this yields (modulo a constant term) the Hamiltonian of the isotropic Landau-Lifshitz ferromagnet With the help of the Cayley-Hamilton polynomial for 3 × 3 matrices, which can be readily integrated The first two terms in this expression exactly match the form of the N = 2 case above, apart from a different multiplicative factor in front. Somewhat unexpectedly, the symplectic generator involves an imaginary term which only reduces to a real quantity in the τ → 0 limit.
Case N = 4. We conclude our analysis by considering also the N = 4 case, which represents the first instance where the generator involves two functionally independent matrix invariants, namely s 2 and s 4 . Let us specialize to the traceless case with s 2j−1 = 0, i.e. signature Σ (2,4) . Two independent equations in (B.26) are 8 ∂H Invoking once again the Cayley-Hamilton theorem, this time using that the odd trace invariants are all zero, we arrive at the recursion of the form The solution can still be found in closed form, whereas s 2m+1 = 0. Performing the resummation, we find a coupled system of two PDEs where F is a constant invertible complex matrix and Ad A (B) = A B A −1 .
Semi-discrete limit. In the continuous time limit, τ → 0, the symplectic maps reduces to the following equation where the Poisson bracket is defined as 3) takes the role of a compatibility condition for an auxiliary linear problem in the form of a semi-discrete zero-curvature condition, see (2.98) in the main text.
In the derivation below we assume the absence of twist, F = 1 (equiv. B = 0) which can be easily incorporated back at the very end. The first step is to promote the Poission structure (C.5) to the level of Lax matrices which yields the quadratic Sklyanin bracket where we have used a short-hand As a consequence of M 2 = 1, we have the following swap operations Exploiting the above identities we have which readily yields a (non-standard) semi-discrete zero-curvature representation with 'left' and 'right' temporal propagators The last step is to transform Eq. (C.14) into the standard form and determine the temporal component of the Lax pair V , +1 (λ). This can be accomplished with aid of the following 'inversion identities' for Lax matrices which combined imply the 'exchange relation', Using these, the right propagators appearing in L (λ)V R +1 (λ) − V R (λ)L (λ) in Eq. (C.14) can be brought to the left in the following manner: The first two terms of the above expression are already in the desired form. The sum of the two Lax operators in the last two terms is composed of a term proportional to the identity and a term that is the inverse of A −1 and will cancel when the last two terms are summed. It is therefore only the term proportional to the identity that will contribute. By inserting the split identity, 1 = −(λ 2 + 1) −1 L (−λ)L (λ), the last expression can be rewritten as By substituting this result in the non-canonical zero-curvature condition (C.14), we eventually restore the standard form (C.17), with the temporal component V (λ) reading where we have simultaneously reinstated the magnetic field.
Zero-curvature formulation in continuous space-time. Having derived the semi-discrete version of the zero-curvature representation, we are in a position to take the continuum theory limit and obtain also the zero-curvature representation of the PDE obeying the continuous version of the zero-curvature condition

D Local conservation laws
Isospectrality (see Section 2.4.1) is one of the central implications of the Lax (zero-curvature) property, resulting in O(L) functionally independent conserved phase-space functions in involution. In integrable systems with local interaction one can typically find local conservation laws, i.e. conserved quantities that can be represented as a spatially homogeneous sum of densities with a compact support. The most common procedure is to expand the logarithm of the transfer map as a power series in λ around a distinguished point λ 0 . In the context of integrable matrix models introduced in Section 2, natural candidates for such expansion points are where one of the Lax matrices L(λ) or L(λ + τ ) in the staggered monodromy M(λ, µ; {M }) degenerate into a projector, with a rank-k projector P (+) ≡ P and its orthogonal complement In what follows we assume k ≤ N/2, in which case the Lax matrices degenerate into P (+) . Employing the projector realization (2.18) with Eq. (2.27) and evaluating the transfer matrices at the appropriate projection points, we obtain where P 0 and g are defined with Eq. (2.18). The logarithm of these expressions however does not reveal any organizing principle which would manifestly generate strict locality. The exception are only the k = 1 cases where the transfer matrices evaluated at λ 0 = −i yields completely factorizable scalar expressions which are considered below.
Local conservations laws in projective models. Here we assume k = 1, that is integrable matrix models with complex projective spaces CP N −1 as target spaces. The associated Lax matrices at two distinguished points λ 0 = ±i degenerate into rank-1 projectors Let us remind that spatial subscripts, such as in Ψ ± , designate that the quantity depends on the local matrix variable M . As a consequence, the transfer maps completely factorize into sequence of scalars (matrix elements) (D.7) In the above expressions the spatial index is understood modulo L, whereas the nomenclature 'even' and 'odd' here refers to the sub-lattices at which the degeneracy occurs.
As an immediate corollary of this construction, the logarithms of the square moduli of the transfer matrices yield conserved quantities which are manifestly local and real, An infinite tower of higher local conservation laws, labelled by integer n ∈ N, can then be produced in an iterative fashion by means of logarithmic differentiation, (D. 15) In this way, we obtain two inequivalent sequences of local charges with densities supported on 2n + 1 adjacent lattice sites. 16 It is worthwhile stressing here that there is no connection between the lowest local conservation law and the symplectic two-body propagator. Indeed, the density of lowest local conservation law is supported on three adjacent lattice sites.
Hamiltonian limit. In the continuous time limit τ → 0 the degeneracy points λ 0 = {−i, −i − τ } collide with one another and the two towers of local conservation laws merge together. This means in effect that an infinite family of local conservation laws can now be generated from the logarithm of the homogeneous transfer map and λ-derivatives thereof. In particular, for the lowest-degree conservation laws we find, (D.17) 16 The outlined construction is only meaningful as long as the support of the densities does not exceed the length of the system, i.e. before 'wrapping effects' take place.
(ψ 0 , ψ 1 , . . . , ψ N −1 ) T generated by some HamiltonianĤ. If we regard |Ψ(t) as a variational wavefunction, its time-evolution corresponds to extremizing a classical action S = dt L[Ψ] with Lagrangian L[Ψ(t)] = Ψ(t)| i∂ t −Ĥ |Ψ(t) . This approach is called the time-dependent variational principle. In practice it is convenient to operate with an unnormalized wavefunction and treat the complex-conjugate field componentsψ j as independent variational variables. With this in mind, we consider a Lagrangian of the form To extract the semi-classical evolution generated byĤ we use generalized coherent states. For simplicity we restrict our considerations first to complex projective spaces CP N −1 = SU (N )/S(U (N − 1) × U (1)) where, in close analogy to the well-known spin-coherent states of CP 1 ∼ = S 2 , we can define the CP N −1 coherent states by the 'vacuum rotation' where g ∈ G = SU (N ) and |0 is the highest-weight state of an irreducible finite-dimensional representation of g = su(N ). We shall consider here only the fundamental representation. While a general group element g is fully determined by N 2 − 1 real parameters (e.g. Euler angles), the number of independent parameters which parametrize |Ψ is actually smaller. This stems from the fact that the vacuum state |0 stays intact under SU (N − 1) rotations in the orthogonal complement of |0 . This means that the unit complex vector |Ψ lies on the real sphere S 2N −1 . However, vectors which only differ by an overall U (1) phase represent the same physical state and must thus be identified, meaning that |Ψ is indeed an element of a coset space S 2N −1 /U (1) ∼ = CP N −1 , a manifold or real dimension 2(N − 1). A general un-normalized variational wavefunction |Ψ ∈ CP N −1 is therefore parametrized by N − 1 complex variables z i , namely local affine coordinates of complex projective spaces, z = (z 1 , z 2 , . . . , z n ) T . The Lagrangian takes the form The later indeed defines a hermitian metric tensor η known as the Fubini-Study metric, and can be generated through differentiation Since the metric tensor is non-degenerate, Det(η) = 0, Eqs. (E.7) can be readily inverteḋ The inverse of the Fubini-Study metric reads explicitly The above construction can be straightforwardly lifted to complex Grassmannian manifolds Gr C (k, N ). In this case the variational wavefunction Ψ becomes a complex matrix of dimension N × k, whereas the Riemann metric tensor Here J > 0 is the ferromagnetic exchange coupling and Π |α ⊗ |β = |β ⊗ |α is the permutation matrix acting in C N ⊗ C N , i.e. to fundamental representations of su(N ). Here and subsequently we shall keep dependence on N implicit throughout the derivation. This class of integrable models, introduced in [85,86], can be diagonalized by means of the Bethe Ansatz [13].
In the basis of traceless hermitian generators X a ∈ g (cf. Eqs. (2.33) and (2.34)), the permutation matrix assumes an expansion where κ ab ≡ 2δ ab .
To obtain the semi-classical limit of Eq. (E.12) we must not allow for any quantum correlations in the variational wavefunction. This amounts to projecting the Hamiltonian onto the subspace of many-body (product) coherent states (E.14) where |Ψ (t) is a CP N −1 coherent state inserted at position . The variational states |Ψ thus belong to the classical phase space M L , pertaining to the low-energy sector of the quantum chain Hilbert space. The task at hand is to derive the quantum-mechanical transition amplitude The first term is a geometric contribution which stems from non-orthogonality of coherent states located at two adjacent time slices, Ψ (i+1) |Ψ (i) = Ψ(t + ∆t)|Ψ(t) = exp − ∆t Ψ(t)|∂ t Ψ(t) + O(∆t 2 ), (E. 19) which, in the ∆t → 0 limit, produces the Wess-Zumino term The kinetic part S kin = dt dx L kin (x) in Eq. (E.18) comes from spin interaction governed by the quantum HamiltonianĤ. To derive the Lagrangian density, let us initial specialize to M 1 ∼ = CP N −1 , where each coherent state can be represented by a rank-1 projector P = |Ψ Ψ|. For computational purposes we instead define a hermitian traceless matrix Y = N g |0 0| g † − 1 N = N P − 1, (E. 21) which is subjected to the non-linear constraint The expectation value ofĤ in a many-body semi-classical state (E.14) is calculated as using the coherent-state averages of the su(N ) generators X a ≡ Ψ | X a |Ψ = Tr(X a P ) = 1 N Tr(X a Y ). (E. 24) In order to take the continuum limit, we first recast Eq. (E.23) in the difference form. Exploiting the completeness relation a,b κ ab Tr(X a A)Tr(X b A) = Tr(A 2 ), valid for an arbitrary traceless N -dimensional matrix A, we have the sum rule 2 a [Tr(X a Y )] 2 = Tr Y 2 = N (N −1).
In this ways we obtain where we have used that 2 a X a X a +1 = 2 a X a 2 − a ( X a − X a +1 ) 2 . With the aid of the long-wavelength expansion (assuming the lattice spacing ∆ = 1/L) we can now pass to the field-theory regime. At the leading order O(∆ 2 ) we deduce using short-hand notion Y x ≡ ∂ x Y (x) and similarly for higher partial derivatives. Finally, replacing the sum with an integral dx and introducing the renormalized coupling J cl = J∆ 2 /2, we arrive at a neat final result Ĥ = J cl N 2 dx Tr Y 2 x = J cl dx Tr P 2 x . (E.28) The kinematic contribution to the Lagrangian density is therefore L kin (x) = J cl Tr((∂ x P (x)) 2 ). To give a simple example, in the case CP 1 ∼ = S 2 we have Y (x) ≡ M (x) = σ · S(x), and using 1 2 Tr(M 2 x ) = S x ·S x we retrieve the Hamiltonian of the isotropic Landau-Lifshitz magnet, with the equation of motion Even thought an explicit construction of an SU (N ) coherent state restricted to Gr C (k, N ) is slightly more involved, the expectation values of the Lie algebra generators still corresponds to the moment maps, that is traces with respect to rank-k projectors P . Up to a shift and rescaling, the latter is equivalent to the traceless hermitian matrix variable Y = N P − k 1 N , obeying Y 2 = (N − 2k)Y + (N k − k 2 )1 N . Repeating the above derivation we once again arrive at Eq. (E.28), apart from a constant shift.