Superuniversality of superdiffusion

Anomalous finite-temperature transport has recently been observed in numerical studies of various integrable models in one dimension; these models share the feature of being invariant under a continuous non-abelian global symmetry. This work offers a comprehensive group-theoretic account of this elusive phenomenon. For an integrable quantum model invariant under a global non-abelian simple Lie group $G$, we find that finite-temperature transport of Noether charges associated with symmetry $G$ in thermal states that are invariant under $G$ is universally superdiffusive and characterized by dynamical exponent $z = 3/2$. This conclusion holds regardless of the Lie algebra symmetry, local degrees of freedom (on-site representations), Lorentz invariance, or particular realization of microscopic interactions: we accordingly dub it as superuniversal. The anomalous transport behavior is attributed to long-lived giant quasiparticles dressed by thermal fluctuations. We provide an algebraic viewpoint on the corresponding dressing transformation and elucidate formal connections to fusion identities amongst the quantum-group characters. We identify giant quasiparticles with nonlinear soliton modes of classical field theories that describe low-energy excitations above ferromagnetic vacua. Our analysis of these field theories also provides a complete classification of the low-energy (i.e., Goldstone-mode) spectra of quantum isotropic ferromagnetic chains.

A complete characterization and classification of dynamical properties of isolated interacting quantum manybody systems remains one of the central unsettled problems in statistical mechanics, whether classical or quantum. Especially in one-dimensional systems, a range of exotic dynamical phenomena have been demonstrated, both theoretically and experimentally. Two prominent examples are integrable and many-body localized systems [1][2][3], which feature extensively many conserved quantities and therefore can persist in nonthermal "generalized Gibbs states" that are measurably different from the standard thermal ensemble [4][5][6][7][8][9]. Because these extensive conservation laws lead to nonstandard equilibrium states, and because hydrodynamics begins with an assumption of local thermal equilibrium, it follows that hydrodynamics is also modified for integrable systems. Thus, instead of normal diffusion, non-disordered integrable systems typically exhibit ballistic transport with finite Drude weights [10,11], whereas in localized models transport is entirely absent [2].
Integrable systems feature coherent quasiparticle excitations with infinite lifetime [12,13] that propagate through the system in a ballistic manner while scattering elastically off one another. The same picture remains valid in thermal ensembles at finite temperature, where one can think of quasiparticles being 'dressed' due to interactions with a macroscopic thermal environment [14]. Thermal fluctuations are responsible for screening, and thus conserved charges carried by quasiparticles can be appreciably different from the bare values. This effect is captured by the versatile framework of generalized hydrodynamics (GHD) [15,16]. Among other results, GHD has led to the explicit characterization of ballistic transport [17][18][19][20] and analytic treatments of various other nonequilibrium protocols [21][22][23][24][25]. Remarkably, despite the ballistic motion of individual excitations, certain integrable models do not necessarily exhibit ballistic transport of macroscopic scales, but instead display normal diffusion or even anomalous diffusion; this is the case for a distinguished subset of conserved quantities linked with internal degrees of freedom, whereas other conserved quantities (such as energy) undergo ballistic transport [26].
This work is dedicated to anomalous transport in integrable models. This unexpected phenomenon was first found in Ref. [27] in the Heisenberg spin-1/2 chain, which can be regarded as an archetypal example of a quantum many-body interacting system. Nowadays there is numerical evidence [28] that the dynamical exponent and asymptotic scaling profiles of dynamical structure factors belong to the Kardar-Parisi-Zhang (KPZ) universality class [29]. Specifically, the diagonal dynamical spin correlations C jj (x, t) ≡ S j (x, t)S j (0, 0) among the spin components S j (for all j ∈ {x, y, z}), evaluated in thermal equilibrium at half filling, have been found to comply with the scaling form characterized by scaling function f KPZ (tabulated in [30]), KPZ nonlinearity coupling parameter λ, and normalized by diagonal static charge susceptibilities C jj = dx C jj (x, t). Kardar-Parisi-Zhang physics is a widespread phenomenon in stochastic growth processes and dynamical interfaces [31][32][33][34]. Its occurrence in deterministic Hamiltonian dynamical systems is therefore an interesting curiosity. The microscopic mechanism responsible for superdiffusion is not yet understood in full detail. The fact that spin (or charge) superdiffusion of the KPZ type has also been numerically observed in a number of other quantum chains, such as the higher-spin SU(2) integrable chains and the SO(5)-symmetric spin ladder [35], the Fermi-Hubbard model [36], and even in the classical Landau-Lifshitz equation [37][38][39] (and its higher-rank analogues that exhibit symmetry of nonabelian unitary groups [40]), give a strong indication that there is a general principle behind superdiffusion in integrable systems with non-abelian symmetries that awaits to be uncovered.
Generalized hydrodynamics (GHD) has already provided a number of invaluable theoretical insights into this question. Specifically for the case of the quantum Heisenberg spin chains, both the dynamical exponent z = 3/2 [41] and the KPZ coupling constant λ KPZ [42] (though not the scaling function) have been inferred with aid of a heuristic extension of the GHD framework, using full advantage of the exact knowledge of the Bethe ansatz quasiparticles. Recently, a phenomenological explanation for the observed KPZ phenomenon has been given in [43] which invokes the notion of hydrodynamic 'soft modes' coupled to an effective noisy environment. The current understanding is that such soft models are a manifestation of the so-called giant quasiparticles in the long-wavelength regime, whose emergent dynamics is governed by a classical action [42]. This appears to suggest that that the observed superdiffusive spin dynamics in classical rotationally symmetric (e.g. SO(3)invariant) spin chains has the same microscopic origin as that of quantum spin chains. In spite of this theoretical progress however, a fundamental question nonetheless remains unsettled: In integrable Hamiltonian dynamical systems, what are the necessary and sufficient conditions for the superdiffusive charge dynamics characterized by anomalous dynamical exponent z = 3/2?
We dedicate this paper to answering this question on general grounds, both in the realm of integrable lattice models and integrable quantum field theories.
Before delving into mostly technical aspects, we would like to offer a broader perspective on the problem at hand. By invoking the universal GHD formulae for the spin diffusion constants, it is evident from the outset that divergent charge diffusion constants are only possible in systems with infinitely many quasiparticle species. This necessary condition is quite generally fulfilled in integrable lattice models, including in particular the models with non-abelian symmetries that will be our primary focus. Given that such models do have infinitely many quasiparticle species, however, there are two scenarios that seem particularly plausible: (I) All integrable Hamiltonians that are symmetric under a non-abelian Lie group G display universal superdiffusive dynamics of the Noether charges in G-invariant (i.e. unpolarized) Gibbs states, (II) integrable models accommodate for a wider range of dynamical exponents, depending possibly on the rank of type of Lie algebras and their representations assigned to local degrees of freedom.
There are several recent studies, e.g. [35,36,44], which speak in favor of (I). Another strong piece of evidence comes from a recent study [40] of classical integrable matrix models, providing integrable spacetime discretizations of higher-rank analogues of SU(n)symmetric Landau-Lifshitz field theories on complex projective spaces and Grassmannian manifolds. In Ref. [40] the authors provide clear numerical evidence for the KPZ scaling profiles independently of rank r = n − 1 and further conjectured that the phenomenon occurs across all classical non-relativistic G-invariant integrable field theories with hermitian symmetric spaces as their target manifolds. Indeed, the class of models considered in Ref. [40] govern low-energy spectra of integrable SU(n)-invariant Lai-Sutherland ferromagnetic spin chains [45,46] that are included as a part of this study. On this basis, it is reasonable to anticipate that coherent semiclassical modes analogous to the aforementioned giant quasiparticles to emerge as a general feature of integrable quantum chains. Their microscopic description is nevertheless not known at this moment.
Nonetheless, option (II) cannot be a priori rejected either. We note that within the phenomenological framework of nonlinear fluctuating hydrodynamics (NLFHD) [47,48], an infinite family of anomalous dynamical exponents can arise [49]. Now, NLFHD does not directly apply to integrable systems, and there are good reasons to doubt that it can capture superdiffusion in the integrable case [43]: notably, in NLFHD, superdiffusion arises as a correction to ballistic transport, rather than as the leading dynamical behavior. Regardless, the existing numerical evidence on dynamical exponents in integrable systems is also ambiguous: for instance, a study of the integrable SU(4) spin chain found a distinct exponent z ≈ 5/3 [50], suggesting the possibility that different nonabelian symmetries may after all realize distinct universality classes of anomalous transport. In addition to its fundamental interest, this question is experimentally relevant because interacting quantum lattice systems possessing SU(N ) symmetry can be implemented using ultracold alkaline-earth atoms [51,52]. In the context of solid state physics, several strongly coupled ladder compounds also realize (approximately) the SU(4)-symmetric ladder in the presence of fields [53][54][55].

Summary
In this work we outline a systematic theoretical analysis of anomalous charge transport, specializing to the class of integrable quantum "spin" chains symmetric under global non-abelian simple Lie groups. In order to provide a universal algebraic description we shall heavily rely on representation theory of quantum groups and the associated fusion relations. For clarity, we here first summarize our main findings.
Our central results is that an anomalous algebraic dynamical exponent z = 3/2 associated with transport of the Noether charges is a common feature of integrable Hamiltonian systems invariant under a general nonabelian Lie group G, provided the equilibrium state does not break the global symmetry by the presence of finite chemical potentials. This statement holds independently of the type of simple Lie algebra and on unitary representations associated with local Hilbert spaces (in quantum chains) or local degrees of freedom (in integrable QFTs): the only requirement is that the charge Q whose correlation functions we study must transform nontrivially under G (unlike, e.g., the energy). We thereby establish superuniversality of superdiffusive charge transport.
Our analysis incorporates the following complementary approaches: (i) We carry out a scaling analysis of the universal Nested Bethe Ansatz dressing equations, concluding that the spectrum of giant quasiparticles (governing the semiclassical long-wavelength dynamics of the charge density in highly excited eigenstates) exhibit the same type of asymptotic scaling relations irrespectively of the non-abelian symmetry algebra. This implies that the kinetic theory argument outlined in Ref. [41] carries through in general. To solidify this conclusion, we derive an explicit closedform solution of the dressing equations for the case of higher-rank unitary groups SU(N ) (including the general dependence on the U(1) chemical potentials coupling to the Cartan charges).
(ii) Secondly, we verify our predictions through tensornetwork based numerical simulations by computing dynamical charge correlations functions for a number of representative cases, including unitary, orthogonal and symplectic Lie groups (shown in Figs. 1 and 2); evidently, all the cases we have studied yield the exponent z = 3/2.
(iii) Lastly, we elucidate the physical nature of the giant quasiparticles. These are none other than semiclassically quantized classical soliton modes, which from the spin chain viewpoint correspond to macroscopically large coherent states made out of interacting, quadratically dispersing (i.e., "magnon-like"), Goldstone modes above a ferromagnetic vacuum. At first glance it may seem counter-intuitive that long-wavelength modes matter to high-temperature physics. Here integrability comes into play, ensuring that these large bound-state excitations remain well-defined quasiparticles even in thermal states; however, their (bare) properties do get dressed by thermal fluctuations. To corroborate this, we establish a one-to-one correspondence between the Bethe ansatz quasiparticles and the spectrum of Goldstone modes. This correspondence is subtle: there are many more distinct Goldstone modes for a given symmetry-breaking pattern than there are magnon species (flavors) in the Bethe ansatz spectrum. Counting the Goldstone modes correctly therefore requires the notion of composite quasiparticles called 'stacks'. To our knowledge, these multiflavored stacks have not been explicitly classified thus far.
While we mostly focus on integrable quantum chains with ferromagnetic exchange, we also shortly discuss in Sec. VI how KPZ superuniversality emerges even in Lorenz-invariant integrable quantum field theories possessing internal isotropic degrees of freedom which take values in compact simple Lie groups G or coset spaces thereof. Such integrable QFTs are characterized by non-diagonal scattering, signifying that their elementary quasiparticle excitations can exchange isotropic degrees of freedom upon elastic collisions. The outcome of that are dynamically produced massless pseudoparticles that are responsible for charge transport at finite temperature. These so-called 'auxiliary quasiparticles' are interacting magnon waves and bound states thereof, which formally resemble magnons of the quantum ferromagnetic spin chains [56]. Divergent charge diffusion constants are thus once again attributed to the presence of interacting giant magnons.
Outline. The paper is structured as follows. In Sec. II we succinctly review the formalism of generalized hydrodynamics, provide closed formulae for the charge diffusion constants, and proceed by summarizing the general physical picture and scaling arguments that lead to superuniversality. In Sec. III we discuss our numerical results on transport. The remainder of the paper consists of a longer technical section IV with all the background information, divided into various subsection devoted to various theoretical aspects, including the details of the Nested Bethe Ansatz diagonalization technique and the algebraic structure of the quasiparticle spectra for a family of quantum spin chains. In Sec. V we elaborate on the semiclassical limit of integrable models and discuss related concepts. In Sec. VII we summarize our results and propose areas for future exploration.

II. OVERVIEW OF METHODS AND RESULTS
We aim to characterize charge dynamics in a Hamiltonian dynamics invariant under a non-abelian continuous symmetry group G. We begin our presentation by first outlining the general setting and introducing the key quantities of the linear transport theory. We shall provide compact spectral representations for the charge diffusion constant in terms of quasiparticle spectra, following largely previous works on the subject [41,[57][58][59][60].
We specialize to simple Lie groups G of rank r, generated by Lie algebra g. Owing to the Noether theorem, the system possess local conserved (Noether) charges, with local densities q (σ) (x), one per each hermitian generator X σ of Lie algebra g. For simplicity of notation we shall not make explicit distinction between lattice and continuum models, i.e. in lattice model spatial integration in Eq. (2) should be understood as a discrete summation. Our first objective is to obtain closed-form expressions for transport coefficients in thermal equilibrium. With no loss of generality we can specialize exclusively to grandcanonical Gibbs ensembles (at inverse temperature β), corresponding to density matrices of the form with normalization (partition function) Z β,h ≡ Tr β,h . Parameters h ≡ {h i ∈ R} are the U(1) chemical potentials which have been assigned to a maximal set of commuting (Cartan) charges Q (i) with i ∈ I r ≡ {1, 2, . . . , r}.
The presence of finite chemical potentials h allows us to study transport at generic values of the charge densities. Their values are of profound importance for the transport phenomena. By adopting generic non-vanishing chemical potentials. When all h i = 0, the ensembles (3) explicitly violate the global symmetry of G, and one is left with a residual symmetry of the maximal abelian subgroup T ≡ U(1) r ⊂ G (called the maximal torus). The complete set of dim(g) Noether charges can be accordingly separated into two sectors: the 'longitudinal' charges Q (i) assigned to the 'unbroken' (Cartan) generators and the non-abelian set of 'transversal' charges Q (σ) for σ / ∈ I r , satisfying [Q (i) , Q (σ / ∈I) ] = 0, associated to the 'broken' generators. In this work we are exclusively interested in emergent anomalous charge transport that arises in the limit h i → 0 where the distinction between longitudinal and transversal correlators disappears (since Ginvariance of β,h gets restored); it will be thus sufficient to focus exclusively to the longitudinal sector (see e.g. [61] for remarks about the nature of transversal correlators).

A. Generalized hydrodynamics
We now formulate the transport theory for integrable systems, known as generalized hydrodynamics [15,16] (cf. [62] for an overview). As our starting point, we begin by an infinite tower of local (including quasilocal [9,63]) conservation laws For the time being k ∈ N is just a formal discrete mode index. Although the complete set of (quasi)local charges can be constructed explicitly with aid of the Algebraic Bethe Ansatz techniques [63], recently adapted in to construct the associated currents [64,65], this step can be in practice circumvented as long as one operates at the level of thermal averages. The GHD formalism provide an explicit prescription to describe large spatio-temporal variations of thermally averaged conservation laws state by Eqs. (4). The key ingredient is to express expectation values of the current densities as functionals of the charge averages. There are various meaningful choices for the thermodynamic state functions one can consider. A particularly useful one is to use quasiparticle rapidity densities ρ A (θ), where the 'type' label A runs over the entire model-specific (thermodynamic) quasiparticle content, and θ is the corresponding rapidity variable parametrizing their bare momenta, that is p = p(θ). An infinite collection of functions {ρ A (θ)} uniquely specifies a macrostate, representing unbiased microcanonical ensembles of locally indistinguishable thermodynamic eigenstates. By exploiting a one-to-one correspondence between the expectation values of the (quasi)local conservation laws and quasiparticle content in an equilibrium macrostate [66,67], it proves useful to recast Eq. (4) as a continuity equation for the quasiparticle densities The quasiparticle current densities take a simple factorized form at the Euler scale [15,16] Here v eff A are state-dependent effective group velocities which determined from the dressed quasiparticles dispersions For a simple Lie group G of rank r, elementary quasiparticle excitations (defined with respect to a reference ferromagnetic order parameter) come in exactly r different types. We shall call these 'flavors'. In addition to that, quasiparticles participate in formation of bound states, a quasiparticle that carries s quanta of flavor a is accordingly assigned a pair of quantum numbers A = (a, s). All of the models we consider have infinitely many species of bound states, i.e., s takes values in an infinite countable set (typically ranging from 1 to ∞).
Expanding above a reference equilibrium densities ρ A (θ), that is ρ A (θ; x, t) = ρ A (θ) + δρ A (θ; x, t), the hydrodynamic evolution of density (or charge) fluctuations δρ A (θ; x, t) on Euler scale is encoded in the linear propagator (flux Jacobian) A = ∂j/∂ρ, reading where (for compactness of presentation) we have employed the tensor notation by flattening the quasiparticle and rapidity labels. Equation (8) can be diagonalized by performing a basis transformation depending (non-linearly) on Fermi functions n of the underlying equilibrium state, where δφ are the normal modes of GHD (defined uniquely modulo normalization of static covariances). In the basis of normal modes, the propagator A acts diagonally, with eigenvalues given by the effective velocities (7), that is The physical interpretation of Ω[n] becomes transparent in the formalism of the Thermodynamic Bethe Ansatz, where it plays the role of a dressing operator for conserved quantities,

B. Diffusion constants
We shall now introduce the linear transport coefficients. Here we define them based on the asymptotic behavior of dynamical structure factors (or dynamical charge susceptibilities). For our intents it will be sufficient to only consider the Cartan charges and introduce the associated r-dimensional matrix of dynamical susceptibilities where bracket • designates the connected correlator evaluated in a grand-canonical Gibbs ensemble given by Eq. (3). Static susceptibilities are accordingly given by the time-invariant sum rule C ij = dx C ij (x, t).
For generic values of the Cartan chemical potentials h i , the variance of the structure factor in integrable models experiences ballistic spreading, signalled by the ballistic dynamical exponent z = 1 and charge Drude weights D ij [10,68,69]. The charge Drude weights admit the following mode resolution [18,19] Heren A (θ) ≡ 1 − n A (θ) denote thermal Fermi occupation functions associated with vacancies (holes). All the thermodynamic quantities in the integrand of Eq. (14) depend on temperature and chemical potentials. Drude weights do not provide complete information about the late-time relaxation of C ij (x, t). To deduce the asymptotic behavior of C ij (x, t) on a sub-ballistic scale, we adopt the kinetic theory framework of Refs. [41,58]. One normally envisions a thermodynamic system divided up into large fluid cells of size , with each cell being approximately in local thermal equilibrium. In a hydrodynamic description, both the dressed charges q (i)dr A and local chemical potentials h i get promoted to dynamical quantities which, in a macroscopic macroscopic fluid cell of length exhibit thermal fluctuations of the order O( −1/2 ), which will, in turn, lead to diffusion.
In generic local equilibrium states, diffusion is a subleading correction to the ballistic transport characterized by the Drude weight. However, in unpolarized thermal ensembles in systems with non-abelian symmetries, the charge Drude weight vanishes and the leading transport behavior is sub-ballistic. Charge Drude weights are proportional to dressed charges q (i)dr A carried by quasiparticles, cf. Eq. (14). The latter are quite different from their bare (quantized) values q (i) A and depend non-trivially on chemical potentials of the background equilibrium state, including the U(1) chemical potentials h i . In unpolarized thermal states that exhibit full invariance under G, the dressed quasiparticle charges vanish simply by symmetry under G (see, e.g., Refs. [17,19,43]): one can see this as the screening of the quasiparticle charge by the thermal environment. Therefore, the charge Drude weight (14) vanishes.
The leading response therefore occurs at the diffusive scale, where one treats the chemical potentials as dynamically fluctuating quantities, with fluctuations that are suppressed by the hydrodynamic scale, O( −1/2 ). For sufficiently small h i the quasiparticles carry dressed charges linearly proportional to h i , i.e., they behave paramagnetically. Fluctuations of chemical potentials induce fluctuations of thermally dressed charges in accordance with Notice that that chemical potentials can be simply related to densities of the Cartan charges via Here C ij = ∂q (i) /∂h j = (∂ 2 /∂h i ∂h j ) log Z are static charge susceptibilities, and it will be helpful below to express them in terms of a mode expansion [18,19] analogous to Eq. (14): The latter also determine the magnitude of charge fluctuations, q (i) q (j) = C ij / . By combining these two results, thermal fluctuations of dressed charges (or dressed susceptibilities) carried by individual quasiparticle modes can be expressed in the form where Physically, Υ ij A can be interpreted as an effective paramagnetic moment, assumes non-trivial dependence on both quasiparticle quantum numbers A = (a, s).
In diffusive dynamics, the variance of the dynamical structure factors C ij which experiences linear growth at late times, characterized by charge conductivity matrix with static susceptibility matrix C and charge diffusion matrix D [70]. A full expression for the conductivity (Onsager) matrix σ has been derived in [57,59] using the form-factors approach, and its diagonal part in [58] using a kinetic theory formulation. Here we provide a compact expression for D in models of higher-rank symmetry (restricted to the Cartan sector), specializing to the limit of vanishing chemical potentials, h → 0. To this end, we substitute the fluctuation relation (18) into Eq. (13) with Eq. (14) on a characteristic scale set by quasiparticles' effective velocities, = |v eff A (θ)|t [41], yielding the following spectral resolution of the conductivity matrix The diffusion matrix is thus proportional to the identity, corresponding to a single value of charge diffusion constant Quantities Υ ij A can be easily computed explicitly at infinite temperature in various integrable spin chains. For instance, in the A n−1 ≡ SU(n) ferromagnetic integrable chains (with onsite degrees of freedom in the fundamental irreducible representation V of dimension dimV = n + 1), static susceptibility matrix C An evaluated at h → 0, written in the non-orthonormal basis H i of Cartan subalgebra t (with Killing metric Tr(H i H j ) = (κ An ) ij , cf. appendix A for details and conventions), reads C An = κ An /(n + 1), and with aid of prescription (93) we find Υ A1 s = 4 9 (s + 1) 4 , Υ A2 a,s = 1 16 (1 + s) 2 (2 + s) 2 , while for n ≥ 4 it depends on the flavor label a ∈ {1, 2, . . . , n − 1}. Its precise form will not matter for our purposes, but note the large-s scaling Υ A ∼ s 4 .
By adopting quasiparticle densities as variational objects (see appendix B for the derivation), and comparing the result to Eq. (22), one finds the following relation between the dressed charge fluctuations and the dressed differential scattering phases K dr AA , We remark that this expression can thought of as a generalization of the so-called 'magic formula' of Ref. [60], which was originally introduced as a way to relate the different expressions of the diffusion constant of the XXZ spin chain with ∆ > 1 obtained from form factors [59], and from a kinetic theory argument [41] similar to the one we used above.

C. Emergence of superdiffusion
From the considerations above, one might expect to find normal diffusion with z = 2, with diffusion constant given by Eq. (22). Explicit computation shows, however, that this result diverges as h i → 0 [41,44]. This divergence has been catalogued for various specific cases; here we will explain why it holds in general. To this end we define the regularized diffusion constant by truncating the spectral sum (22) at some string index s * , where D a,s ≡ D A corresponds to each summand in Eq. (22). As a finite sum, Eq. (26) is manifestly finite. If D * converges as s * → ∞ then one has regular diffusion. Suppose instead (i) that D * diverges as D * ∼ (s * ) υ . This already signals superdiffusive transport. We further posit (ii) that the time-dependent diffusion constant D(t) ∼ t α with 0 < α ≤ 1 and (iii) that the dressed velocities of quasiparticles fall off at large s as v eff s ∼ s −ν with exponent ν > 0. We now derive scaling relations among these exponents. We note first that the distance over which a spin packet has spread in time t is defined by x a (t) ∼ D(t)t ∼ t (1+α)/2 , thus the dynamical exponent z ≡ 2/(1+α). At time t we can divide quasiparticles into "light" quasiparticles for which v s t ∼ t/s ν > x a (t) and "heavy" quasiparticles for which v s t < x a (t). The natural crossover scale at time t is thus set by t/s ν * ∼ t (1+α)/2 , so s * ∼ t (1−α)/2ν . Combining this with assumption (i) we conclude that In the isotropic Heisenberg chain, analyzed in ref. [41], one has specifically υ = ν = 1 and hence z = 3/2. This type of reasoning appears to suggest that more generally a larger set of anomalous dynamical exponents could be realized in integrable models by appropriate tuning of ν and υ. However, ν and υ are not truly independent and unrestricted parameters in integrable models, since υ depends on the scaling behavior of Eq. (26), which in turn depends on ν through Eq. (22).
This imposes quite stringent and general constraints on υ, as follows. Notice that as h → 0 the static charge susceptibilities C ij , given by Eq. (17), must approach a constant value. At finite h, the occupation factor for a quasiparticle specie a of type s (a bound state of s magnons) is suppressed with a factor of exp(−s h a ), and thus quasiparticles with s h a 1 do not contribute to the susceptibility at chemical potential h a . The label 'a' plays no role in the following argument, and so we will drop it. Writing uniformly h a → h, truncating the sum (17) at s ∼ 1/h, and using the paramagnetic behavior of the dressed charges, we have C ij ∼ h 2 r a=1 s<1/h C ij a,s . Requiring that the C ij have a nonvanishing limit as h → 0, we see that the sum over s has to scale as ∼ s 2 , so the summand C ij a,s must grow linearly in s. The summands in Eqs. (17) and (22) are the same, except for a factor of |v s | ∼ s −ν . Therefore, the diffusion constant must scale as D * ∼ s≤s * s 1−ν ∼ s 2−ν , yielding the scaling relation Finally, the analytic structure of Bethe ansatz solutions implies that exponent ν ≥ 0 is integer-valued. The three remaining possibilities are ballistic transport (ν = 0), which is ruled out by the vanishing Drude weight discussed above; KPZ scaling (ν = 1), which we will argue is generic; and diffusion with possible logarithmic corrections (ν ≥ 2), which occurs in certain fine-tuned models [71] where the single-magnon dispersion is fine-tuned to scale as ω ∼ k 3 (or slower) at long wavelengths k → 0.
In what follows we will establish z = 3/2 within a universal algebraic description of the thermodynamic dressing equations and link them to the underlying symmetries structures and representation theory of quantum groups (Yangians). We find, remarkably, that the Fermi functions assume universal algebraic scaling at large-s, which comes out as a direct corollary of fusion identities amongst the quantum character associated to the Yangian symmetry [72][73][74][75][76]. Similarly, we find that the total state densities and the dressed velocities of giant magnons (when multiplied by regular function and integrate over the rapidity domain, cf. Sec. IV G) decay as respectively, for all flavors a = 1, . . . , r. Most remarkably, these large-s scaling properties hold irrespectively of the local on-site degrees of freedom (i.e. finite-dimensional irreducible unitary representations of g); they can be understood as kinematic constraints stemming from the underlying quantum symmetry algebra. Finally, notice that Υ A ∼ s 4 in conjunction with the above scaling relations implies υ = 1. The upshot is thus that the anomalous fractional algebraic dynamical exponent z = 3/2 is deeply rooted in the fusion relations for the quantum characters.

III. NUMERICAL ANALYSIS
We verify our predictions for the superdiffusive charge transport on a number of representative instances of integrable quantum spin chains invariant under simple Lie groups from the classical series. Specifically, we consider homogeneous spin-chain Hamiltonians in terms of fundamental degrees of freedom given by Eqs. (45). Exceptional algebras are not included in our analysis. We moreover leave out the B-series as the fundamental integrable B 2 ≡ SO(5) chain has been already been studied numerically in a recent paper [35].
We employ numerical tensor-network based computations of the dynamical correlation function C ii (x, t) (as defined in Eq. (1)) in the canonical Gibbs equilibrium at infinite temperature. Owing to the non-abelian symmetry G of both the Hamiltonian H and infinitetemperature Gibbs density matrix, all the components i = 1, 2, . . . , r are identical and we thus suppress the redundant index i. We performed the time-evolving blockdecimation (TEBD) algorithm in the Heisenberg picture [77][78][79], evolving the local charge density operator q(0, t) in a matrix product operator (MPO) representation.
Details. A fourth-order Trotter decomposition is used to propagate the operator forward in time-steps of size δt, followed by truncations at each bond keeping the largest χ Schmidt states. The simulations are accelerated by taking advantage of time-reversal symmetry and translation symmetry to compute the full correlator C(x, t) = q(x, t/2)q(0, −t/2) using a single MPO evolution; additionally, the TEBD implements the maximal abelian subgroup of the on-site symmetry group for efficiency. The computations are checked for convergence in δt and bond dimension χ (up to χ = 1024).
Our TEBD scheme uses a Trotter step size of δt = 0.4. The operators are truncated initially with a constant discarded weight ε = 10 −8 , allowing the bond dimension χ to grow until it reaches a threshold χ max = 512. Subsequent truncations keep only χ max states. The simulations are checked for convergence in the step size δt, the truncation error ε, and the threshold bond dimension χ max . Additionally, the results shown here include rescaling of the correlations at each time to explicitly enforce the charge sum rule dx C(x, t) = C(0, 0), which improves the convergence significantly.
Results. Our main results are succinctly summarized in Fig. 1 where we display the dynamical width of the charge profiles We find clear signature of an asymptotic power-law growth ∆x(t) ∝ t 1/z , and the return probability shows power-law decay C(0, t) ∝ t −1/z , with numerically estimated dynamical exponent z = 3/2 with great accuracy. The scaling collapse in Fig. 2 shows the correlators obey the asymptotic scaling form C(x, t) t −1/z f sc (x t −1/z ). Comparing f sc to the KPZ scaling function f KPZ given by Eq. (1) we find some discernible deviations in the tails (including the basic SU(2) case studied previously in [28]). On accessible timescales, we are unable to infer whether this discrepancy persists at late times or is it merely due to transient effects.

IV. THEORETICAL BACKGROUND
In the rest of the paper we provide the necessary technical details using the language and framework of the Algebraic Bethe Ansatz. Our aim is mostly to give a concise description the quasiparticle content and a grouptheoretical formulation of the TBA equations. This will require us to briefly revisit the key notions from the theory quantum integrability. Despite a vast body of work on the subject, we are not aware of a self-contained exposition of the Nested Bethe Ansatz (NBA) techniques and its various mathematical underpinnings; even some of the fundamental results appear to be dispersed across several specialized articles. We would thus like to use this opportunity to partly fill this gap and offer a comprehensive exposition aiming at physicists with interest in the field of many-body statistical physics out of equilibrium. Our analysis crucial relies on several basic notions from the representation theory of quantum groups known in the literature as Yangians. Specifically, we shall rely on the knowledge of functional identities amongst their 'quantum characters'. The study of Yangian modules has been initiated by Kirillov and Resthetikhin in their seminal papers [80][81][82]. Many key developments in this field are summarized in the review article by Kuniba et al. [83] which will be extensively used.
Concerning physics applications, it is worthwhile mentioning that quantum character formulae proved invaluable for deriving the dilogaritm identities for the conformal central charges (that can be extracted from the lowtemperature scaling limit for a wide range of integrable quantum chains), see e.g. Refs. [83][84][85][86]. We are nonetheless not aware of any applications in the nonequilibrium physics department. Readers not interested in the details of the Bethe ansatz solutions to these models can feel free to skip this more technical section, although some of the scaling properties we derive in subsections IV F and IV G are central to our analysis; the results derived in these two subsections are, to our knowledge, new.
For the most part we shall concentrate our analysis to quantum lattice models with isotropic interactions, adopting finite-dimensional irreducible representations of non-abelian simple Lie algebras g as local degrees of freedom. These models can be viewed as higher-rank variants Node a and b = a share a bond if K ab = 0. The number of bonds between adjacent nodes equals −K ab , corresponding to different angles between simple roots. For non-simply-laced algebras, the arrow is pointing towards the short root.
of ferromagnetic quantum chains that possess manifest global continuous symmetry of a Lie group G. We shall make extensive use of the group-theoretical language. For a brief reminder on simple Lie algebras we refer the reader to appendix A, where we also fix some of our notations. The complete classification of simple Lie algebras is a landmark achievement due to Cartan. The list comprises four infinite families commonly referred to as the classical series, alongside five exceptional algebras g 2 , f 4 , e 6 , e 7 , e 8 . While these exceptional cases can treated in parallel to the classical algebras, they will be of secondary interest for us. We also remind that there are several exceptional isomorphisms amongst the low-rank algebras: Irreducible representations. Finite-dimensional irreducible representations V Λ of g are specified by r nonnegative integer labels m a , called Dynkin labels. They are coordinates in the 'ω-basis' of the fundamental weight ω a . We employ notation Λ = [m 1 , m 2 , . . . , m r ], or sometimes write simply (n), referring to the dimension of V Λ . Writing m i inside the corresponding node of the Dynkin diagram we obtain a one-to-one correspondence between the enumerated Dynkin diagrams and the finitedimensional irreducible representations of V Λ , as shown in Fig. 4. More more information we refer the reader to appendix A.
To every irreducible representation of the highest weight we can also bijectively associate a Young diagram (partition). The latter consists of r rows, with each a boxes in each row such that moving from top to bottom the number of boxes never increases. The upper row contains 1 = r i=1 m i boxes, while in every subsequent row the number is lowered by successively subtracting Dynkin labels, namely a = a i=1 m i . This means that the ath fundamental weight ω a has a single non-zero Dynkin label 1 in the ath node. The associated diagram is a single column with a boxes. It is important to stress that rectangular irreducible representations play a distinguished role as the constitute a closed set of fusion relations. In the context of quantum integrable models solvable by the Bethe Ansatz, rectangular diagrams indeed correspond to distinct quasiparticles species that can be excited in the thermodynamic eigenstates.

A. Integrable lattice ferromagnets
We now proceed with the algebraic foundations of the Bethe Ansatz. We shall describe two complementary principles which underlie quantum integrability. From a formal perspective, integrability arises from solutions to the quantum Yang-Baxter equation (YBE) [87] Here matrices R ij = R ij (u i , u j ) act non-trivially on two copies of vector spaces V i and V j (irreducible g-modules of Lie algebra g) and depend additionally on a pair complex parameters u i,j ∈ C. Geometrically speaking, equation (33) therefore establishes the equivalence of two a priori distinct intertwining protocols on a three-fold space V 1 ⊗V 2 ⊗V 3 . In algebraic terms, R-matrices provide 'structure constants' for the defining relations of quasitriangular Hopf algebras, widely known in the literature as quantum groups [72][73][74]88]. For every simple Lie algebra G, the YBE admits a class which are rational functions of the difference parameter u = u i −u j . The associated quantum algebras are known as Yangians (see e.g. Refs. [72,73,88]).
Equation (33) has a suggestive physical meaning when viewed as an equivalence relation for a factorizable threebody scattering processes, signifying that two apriori different three-body scattering events decompose in terms of sequential two-body scatterings in two equivalent ways. This powerful symmetry principle is intimately tied to presence of coherent long-lived interacting quasi-particles that undergo completely elastic scattering [13]. In this respect, the variables u i parametrize quasi-particle momenta.

B. Commuting transfer matrices
Integrability most prominently manifests itself through the existence of infinitely many commuting operators conserved under time evolution. Let H = V ⊗L Λp denote the Hilbert space of a homogeneous quantum chain of length L, with an irreducible onsite representation V Λp . The row transfer matrices T Λa (u) = Tr Λa M Λa (u) on H, i.e. auxiliary traces of monodromies provide the generating operators for an infinite family commuting charges. Above we have used a standard notation R Λa,i denote the embedded R-matrices which operate non-trivially only on the auxiliary space V Λa and physical lattice site i. Mutual commutativity for u, u ∈ C and for any two arbitrary finite-dimensional gmodules of the highest weight type Λ a and Λ a , namely follows as a direct corollary of (leaving dependence on u implicit) The latter lifts the local property (33) to the entire manybody Hilbert space H. The simplest R-matrices operates on two copies of fundamental g-modules corresponding to the one-box tableaux, V ω1 ≡ V , and we accordingly denote them by R , (u).
Other R-matrices, acting in higherdimensional (physical or auxiliary) representations, can be systematically constructed via an appropriate fusion procedure. This procedure is described in detail in e.g. [89,90]. Amongst T Λa (u), we find an infinite subset of transfer matrices T a,s (u) with rectangular auxiliary representations of the highest weight Λ s ωa ≡ Λ a,s , with with Dynkin labels Λ a,s ≡ [s, s, . . . , s] (depicted by a×s Young diagrams). From an algebraic perspective, they are distinguished by the property that they constitute a closed set of functional fusion identities in the form of Hirota relations. On the other hand, their representation labels (a, s), with 1 ≤ a ≤ r and s ∈ N, bijectively enumerate the types of quasi-particle excitations in the their thermodynamic eigenstates.

Quantum R-matrices
In this section we outline the construction of commuting transfer matrices for homogeneous quantum chains with degrees of freedom in the fundamental representation. We shall thus only operate with the fundamental R-matrix R , (u). For compactness of notation, we subsequently denote simply by R(u).
We begin by SU(n)-symmetric R-matrices which have originally appeared in papers of Lai [45] and Sutherland [46]. The fundamental ones, acting on the product of V ∼ = C n , read simply where denote permutation operators acting in C n ⊗ C n , with unit (Weyl) matrices (E ij ) kl = δ ik δ jl . The orthogonal and symplecic cousins of (37) have been given by Resthetikhin [91]. They contain an extra non-invertible element Ξ = Π T1 = Π T2 , that satisfies In particular, the O(n)-invariant R-matrix acting on two fundamental (vector) representations has the structure Finally, in the case of symplectic groups Sp(2r) (r ≥ 2), the Sp(2r)-symmetric fundamental R-matrix acting in the two-fold product space C 2r ⊗ C 2r reads satisfying All the above R-matrices fulfil the unitary condition R(u)R(−u) = γ(u)1, for some appropriate scalar function γ(u).

Spin chain Hamiltonians
An infinite tower of mutually commuting Hamiltonians with k-site local densities h (k) can be produced by differentiation of the logarithm of T (u), evaluated at a special point u = 0 (where R(0) Π). The (k − 1)-st Hamiltonian in the family has a two-body density h (k) , corresponding to some G-symmetric interaction involving k adjacent fundamental degrees of freedom (with Λ p = ), for k ≥ 1. We subsequently consider the Hamiltonians H ≡ H (2) with nearest neighbor interactions h ≡ h (2) . Depending on g, the Hamiltonian densities h assume the following form

C. Nested Bethe Ansatz
Completely elastic scattering. The algebraic principles of quantum integrability can be neatly translated into the scattering theory language. The many-body Smatrix, albeit being non-trivial, fulfils extremely stringent constraints owing to infinitely many local conservation laws; they ensure only the usual total conservation of particle momenta, but also that the sets of incoming and outgoing momenta are identical. In other words, quasiparticles every pair-wise scattering retain their momenta and consequently the entire many-body scattering matrix completely factorizes into a sequence of two-particle scattering events.
To demonstrate the basic principle, let us consider a one-dimensional system system enclosed in a finite compact region of space of circumference L. Requiring periodicity of the wave-function yields a simple constraint providing the quantization condition for quasi-particle momenta p j . Since scattering is purely elastic and there is no particle production or decay, we only had to take into account that upon every collision a quasi-particle experiences a phase shift encoded in a momentumdependent U(1) scattering amplitude S(p j , p k ). Equation (49) evidently resembles that of free quasiparticles, 'corrected' by the net effect of interparticle interactions.
For the class of model under consideration, we deal with a more general situation where quasi-particles carry additional quantum numbers coming from isotropic degrees of freedom. Condition (49) has to be accordingly promoted to a matrix equation, where now S denotes a scattering matrix, and |Ψ a many-body wave-function. By diagonalizing this resulting equation, the quantization condition takes the form of a coupled scalar equations known as the nested Bethe Ansatz (NBA) equations (cf. [56,[92][93][94]).

Bethe Ansatz diagonalization
The task at hand is to simultaneously diagonalize the set of commuting transfer matrices. In distinction to free systems, many-body eigenstates depend on solutions to the NBA equations (50). There are various routes to tackle this problem, and we shall briefly outline two standard alternative methods for constructing the complete set of Bethe eigenstates.
The best known approach is the celebrated algebraic Bethe Ansatz which mirrors the familiar secondquantization construction. Given a symmetry group G of rank r, one can identify a set of r quasiparticle creation operators B a , one for each simple root of g. The discrete label a can be regarded as a flavor index. The operators B a appear as off-diagonal elements of the fundamental monodromy matrix defined in Eq. (34), see for instance ref. [95]. Identifying a ↔ α a , the Bethe eigenstates assume the form for a suitable set of N a ∈ N complex 'quantum numbers' {θ (a) j } Na j=1 (one for each a ∈ I r ) obtained as distinct solutions to Eq. (50). Coefficients F (a1,...,a N ) indeed turn out to be the wavefunction amplitudes of an integrable quantum chain of rank r − 1 and length N , for which rapidities take the role of inhomogeneities (for further details and a pedagogical introduction the reader is referred to [96,97]). Here we add that one can formulate an alternative algebraic construction of nested Bethe eigenstates by using a single excitation operators for a monodromy matrix twisted by a generic invertible element [95].
The reference state |∅ is a trivial (particle-less) 'Fock vacuum' eigenstate usually called the Bethe pseudovacuum. In the case of homogeneous ferromagnetic chains, it always corresponds to the fully polarized state (i.e. having the maximal weight in H) pointing along a specified direction. Beware that such a state is only unique after gauge-fixing since by virtue of G-invariance Hamiltonians H possess continuously degenerate ground-state manifolds. This signifies that |∅ indeed carries a polarization degree of freedom which takes values on a coset manifold G/H, where the subgroup H ⊂ G is referred to as the stability group (leaves |∅ intact, apart from a phase). With no loss of generality, we have freedom to adopt the conventional gauge choice and set the vacuum polarization to |∅ = L |0 . The Algebraic Bethe Ansatz provides a procedure for explicit construction a complete of the highest-weight Bethe eigenstates |{θ (a) j } , each of which is ascribed a unique set of Bethe roots.
The structure of the eigenvalues of commuting transfer matrices can be most suggestively described in the language of functional Bethe Ansatz, exploiting the fact that transfer matrices can be further decomposed in terms of the so-called Baxter Q-operators (see e.g. [98,99] for further details and explicit construction). Eigenvalues of commuting Q-operators are called Q-functions [98][99][100]. In the nested Bethe Ansatz, eigenvalues of commuting transfer matrices are all functionally dependent on r Baxter Q-functions, given by polynomials of the form Intuitively speaking, the Q-functions can be thought ot as many-body wave-functions in the first-quantization picture whose zeros (called Bethe roots) play the role of 'nodes'. This is made most explicit in the so-called Sklaynin SoV approach [101], where Baxter Q-functions become coefficients in the expansion over basis of (quantum) separated variables [95]. Nested Bethe equations. We consider homogeneous isotropic ferromagentic chains with unitary irreducible g-modules V Λ as the local degree of freedom, specializing to rectangular representations with Λ p = s p ω ap (here subscript p stands for 'physical'). Requiring that a Bethe state Eq. (51) is an eigenstate of the commuting transfer matrices, the cancellation of all the 'unwanted terms' yields a set coupled equations of the form [102] θ (a) wherem a ≡ (Λ p , α a ) = (s p /t a )δ a,ap and (·, ·) is a bilinear form for on the dual space t * of Cartan subalgebra (for definition, cf. appendix A). Equations (53) are none other that the aforementioned NBA equations in disguise. To reconcile them with Eqs. (50), the term in the square brackets on the lefthand side of Eqs. (53) needs to be interpreted as the exponential of the bare momentum, whereas the right-hand side in Eq. (53), expressed abstractly in terms of the Baxter's Q-functions, should be equated with the product of (rational) scattering amplitudes. The physical interpretation of the resulting NBA equations is as follows: spin waves of quasiparticles (magnons) of flavor a + 1 propagate on a (fictitious) lattice formed by quasiparticles of 'adjacent flavor' b whenever there is interaction between the two, namely (α a , α b ) = 0. The NBA equations (53) are in to one-toone correspondence with enumerated Dynkin diagrams, as exemplified in Fig. 4. Fused transfer matrices. The framework of the algebraic NBA allows simultaneous diagonalization of the entire family of fused transfer matrices T a,s ascribed to rectangular irreducible auxiliary representations V a,s . Note that an infinity family of transfer matrices T a,s are not independent objects as they can all be generated out of the fundamental one T via the so-called 'fusion procedure' [89,90]. As a consequence, the infinite family of fused row transfer matrices T a,s (θ) enclose certain functional relations, and (due to their commutativity) likewise for their eigenvalues T a,s (θ). For example, in the su(n) quantum chains, the eigenvalues satisfy a system of coupled recurrence relations of the form in which one can recognize the celebrated Hirota bilinear relation. For compactness of notation we have suppressed dependence of the spectral variable θ and employed a shorthand notation for imaginary half-unit shifts To fix a solution of an infinite hierarchy of T -functions it is thus sufficient to prescribe the 'initial condition' in the form of fundamental eigenvalues T a,1 (θ) which, alongside the trivial boundary functions T 0,s = T n,s = 1, unique determine the remaining T -functions by iterative application of Eqs. (55).
It is worth mentioning that Hirota equations lie at the bedrock of classical and quantum soliton theories [103]. At a formal level, one can think of its as a completely integrable classical lattice gauge theory. The particular version of the Hirota relation has been originally found by Hirota in space-time discretization for a certain class of nonlinear integrable partial differential equations [104][105][106]. Analogous functional relations have been derives for other types of simple Lie algebras g, see e.g. [83] for an extensive review. In the next section we shall outline how the Hirota-type relation reemerges as functional relation associated to certain thermodynamic state functions.
Transfer matrix eigenvalues. The NBA equations (53) are in fact implied by the fusion relations for the commuting fow fused transfer matrices T a,s (θ) associated to rectangular irreducible representations. To outline the logic, we have to switch to the functional Bethe Ansatz perspective. To outline the basic logic, To this end, it is instructive to first shortly examine the nonnested case of the homogeneous SU(2) spin chain with with Q 0 (θ) ≡ θ L . Taking now into account that both T (θ) and Q-functions are manifestly finite-degree polynomials in θ, the roots of Q 1 (θ) must be positioned such to ensure that all the superficial poles vanish. By demanding that residues on the right-hand side of Eq. (56) vanish one recovers precisely the SU(2) Bethe equations.
In the of higher-rank symmetries, eigenvalues T a,s (θ) can be decomposed in terms of r distinct Q-functions. We consider here as an example the case of unitary groups SU(n). By virtue of Eq. (55), it is sufficient to focus on the fundamental transfer matrix T (θ). By formally solving Eqs. (55) in a recursive manner with a chain of Bäcklund transformations, the fundamental eigenvalues T (n) (θ) can be expressed in a closed form [89,107] with boundary functions Q 0 (θ) = θ L and Q n ≡ 0. The NBA equations (50) are implied the pole-cancellation condition, i.e. demanding that T a,s (θ) are entire functions in the θ-plane.
Finite-volume eigenstates. Each highest-weight eigenstate of T (θ) is associated with a degree-L polynomials T (θ). Amongst its L zeros, there arē N 1 real ones (called holes). Provided that Q 1 (θ) has only real roots, one can readily infer from the analytic structure of Eq. (56) for large L (assuming θ ∼ O(L 0 )) using that Q ± 0 Q ∓∓ 1 = T Q 1 for Im(θ) ≷ 0. Complex zeros of T (θ) are approximately roots of Q 1 (θ) shifted in the imaginary direction by ±i. Therefore, when deg Q 1 = N 1 , we can deduce the following equality Ferromagnetic vacuum is 'empty' and fully filled with N 1 = L holes. The opposite of it is the eigenstate with no holes,N 1 = 0 and N 1 = L/2, residing in the maximally filled sector.
A similar analysis applies to the higher-rank models. For example, in the SU(n)-symmetric chains there are n − 1 flavors of magnons and one can now deduce a set of inequalities for a = 1, . . . , r (with N 0 ≡ L and N r+1 ≡ 0). Such 'triangular conditions' are equivalent to the requirement that all Bethe eigenstates are associated a non-negative highest-weight. The maximally saturated sector in the fundamental su(n) chain thus corresponds to filling fractions ν a ≡ N a /L = (n − a)/n. Complex solutions (bound states). The general structure of finite-volume eigenstates is in general far more complicated, mostly because Bethe roots of individual solutions to Eqs. (53) are generically complexed-valued (by reality condition, all the roots must however occur in complex-conjugate pairs). In particular, typical large-volume eigenstates consist of compounds of complex roots with the same real part -known as Bethe strings. The structure of the Bethe equations (53) indicates that for large L (and θ ∼ O(L 0 )), the imaginary separation of two adjacent complex roots within a string-like compound is precisely i, up to finite-size corrections (typically suppressed exponentially in L). Since finite imaginary part in momentum signifies finite localization length, such solutions pertain to multi-magnon bound states. For example, in the above SU (2), an sstring (s ≥ 2) solution induced a complex-conjugate pair of zeros in the corresponding transfer matrix eigenvalue with Im(θ) = ±(s − 1)i.
An analogous analysis can be performed for the higherrank models with nested spectra (see e.g. [108] for a general discussion for the case of unitary superalgebras). Presently, the large-volume Bethe eigenstates in G-invariant ferromagnetic chains can be partitioned in terms of Bethe strings of different flavors a ∈ I r . Specifically, an s-string is made of s constituent magnons with complex rapidities θ (a,s) k modulo finite-size deviations which are negligible in the L → ∞ limit. Notice the exact one-to-one correspondence between the Bethe strings (bound states) and unitary irreducible rectangular representations of g.
For large L, each Q-function Q a (θ) associated with highest-weight Bethe eigenstate can be thus partitioned into strings, , neglecting vanishing finite-size corrections. The total number of sstrings of flavor a, denoted by n a,s , can be computed with an explicit combinatorial formula (given by Eq. (C3) in appendix C 1), counting all admissible rearrangements of 'particles' and 'holes'. We mention in passing that despite such a formula establishes 'combinatorial completeness' of solutions to the NBA equations, it is known the certain legal but irregular solutions involving 'collapsed strings' do occur finite-volume spectrum [109,110]. Such exceptional solutions nonetheless amount to a thermodynamically vanishing fraction and hence are not detrimental to the string hypothesis and the validity of the TBA approach (completeness of the finite-volume NBA equations has been recently proven using an alternative Wronskian formulation in Ref. [111]).
Irreducible multiplets. Owing to G-invariance of the considered Hamiltonians, eigenvectors organize into irreducible subspaces (multiplets) of g. The total number of distinct highest-weight eigenstates equals the number of inequivalent irreducible multiplets in the decomposition of the L-site Hilbert space H ∼ = V ⊗L Λp . The outlined NBA construction only gives access to the highest weight Bethe eigenstates, whereas the descendant eigenstates have to be generated by successively acting with the global lowering generators on every highest-weight Bethe state. In effect, eigenstates within a given multiplet are all characterized by the same set of Bethe roots and consequently yield the same value of (quasi)local conserved charges deriving from T a,s (θ) -they are thus only distinguished by the valued of the Cartan U(1) charges.
Commuting row transfer matrices T a,s (θ) acting on the tensor product Hilbert spaces are not irreducible objects but instead decompose non-trivially into finitedimensional irreducible representations of Yangians Y (g) (commonly known as the Kirillov-Reshetikhin modules [81,82], or W-modules, for short). To be specific, a homogeneous transfer matrix acting in an L-fold tensor product of local Hilbert spaces V Λp corresponds to a reducible W-module W ⊗L Λp . Multiplicities of irreducible components are precisely in agreement with the number of distinct highest-weight Bethe eigenstates. However, when viewed as g-modules, irreducible W-modules generically reducible in a non-trivial manner (with the exception of A n -modules which stay irreducible). Further details on reducibility of Kirillov-Reshetikhin modules are provided in appendix C.

D. Thermodynamic Bethe Ansatz
Having described the structure of finite-volume eigenstate, we proceed to formulate the thermodynamic description of a nested Bethe Ansatz system. This will be achieved in the standard manner by employing the universal functional integral approach introduced originally by Yang and Yang [14], adapting it to accommodate for integrable quantum chains invariant under Lie symmetry group G.
Thermodynamic limit. To infer thermodynamic properties, one has to consider the large-volume limit and retain only eigenstates with an extensive number of excitations, N a ∼ O(L) → ∞, while keeping all the filling fractions N a /L fixed. In doing so, one can afford to omit certain details about individual eigenstates which considerably simplifies analytic treatment. In a generic thermodynamic eigenstate, a typical separation of nearby rapidities (or string centers θ (A) k , in case of bound states) is of the order O(1/L), which makes it natural to introduce an infinite set rapidity densities, one for each thermodynamic quasiparticle specie A = (a, s). Quantities ρ A (θ)dθ representing the total number of quasi-particles of type A residing in an infinitesimal interval [θ, θ + dθ] per unit length. Similarly, one introduces the total densities of available states ρ tot A (θ), encoding densities of the available modes (with mode numbers n i ) in the spectral plane, corresponding to taking the logarithm of the thermodynamic limit of Eq. (53). In contrast with free theories, where the latter takes a constant value of 1/2π, the net contribution of inter-particle interactions induces a non-trivial dependence on a state. It is also convenient to define the densities of holes (i.e. unoccupied modes) . Bethe-Yang equations. Owing to interparticle interactions, quasiparticle densities ρ A (θ) and total densities of available states ρ tot A (θ) are not independent of one another. The relation between the two can be inferred from the thermodynamic limit of the Bethe equations (53), yielding a coupled system of integral (Bethe-Yang) equations [14,112,113] expressing ρ tot A as a functional of ρ A . Here and subsequently, is an abbreviation for the convolution type integrals, which (for dummy functions f A , g A and F AB ) read The two-body kernels in Eq. (62) are given by the logarithmic derivative of the scattering amplitudes S AB (θ), attributed to a two-body scattering event between quasiparticles of types A = (a, s) and B = (a , s ). The complete set of scattering amplitudes S AB can be obtained from the elementary magnon amplitudes via fusion; this amounts to merge a finite 'string' of scattering ampli- Macrostates. Any given admissible set of quasiparticles densities ρ A (θ) (that is a solution to coupled integral equations (62)) is compatible with exponentially many distinct thermodynamic Bethe eigenstates. It is important to stress that individual eigenstates are indistinguishable at the level of local observables [114][115][116], which can be viewed as a manifestation of the 'eigenstate thermalization principle' (in a generalized sense). In this sense, set {ρ A (θ)}, commonly referred to as a macrostate, can be understood as a equivalence classes of thermodynamic eigenstates with a finite density of excitations. In other words, momentum (or rapidity) densities {ρ A (θ)} provide a complete spectral resolution of a local equilibrium state. For the present class of model, every macrostates is ascribed a Fermi-Dirac statistics, i.e. the associated combinatorial weight has the form This expression coincides with the von-Neumann entropy of a generalized Gibbs ensemble associated with a macrostate {ρ A (θ)}.
Equilibrium macrostates can alternatively be (uniquely) parametrized in terms of the expectation values of the local conservation laws [66,67,117]. The caveat here is that when the quasiparticle spectrum involves bound states, the standard set of local conservation laws (obtained from the series expansion of fundamental transfer matrices T (θ)) are actually not sufficient (see Refs. [115,118]) and it is imperative not to leave out the quasilocal conserved charges (generated out of fused row transfer matrices, see Refs. [63,67,117] and a review [9]). For explicit construction and complete classification of generalized Gibbs ensembles exemplified on the emblematic case of the Heisenberg spin-1/2 chain we refer the reader to Refs. [119,120].

TBA equations
We now outline how thermodynamic free energy can be computed in the framework of the Thermodynamic Bethe Ansatz. For definiteness we now confine ourselves to the grand-canonical Gibbs states (3).
Firstly, to every Cartan generator H i ∈ t we associate a conserved Cartan charge Q (i) = dx q (i) (x), with local density q (i) . To fix a basis of the Cartan subalgebra, we define H i ∈ t via H i = [E αi , E −αi ] (for every simple roots α i ), after adopting a particular normalization convention for the Weyl generators, such as the common one κ αi,−αi = 1 (which also fixes κ ij ). The most general element g 0 ∈ T from the torus subgroup can be parametrized as in terms of r scalars x i ∈ R being certain functions of the U(1) chemical potentials h i . Eigenvalues of Q (i) acting on the highest-weight Bethe eigenstates read where q (i) a are (bare) quanta of the ith Cartan charge carried by a quasi-particle of flavor a. Taking further into account that general thermodynamic eigenstates involve bound states (Bethe strings), these further split amongst N a,s Bethe's s-strings (with N a = s s N a,s ) carrying bare charges q (i) a,s . The defining property of macrostate densities {ρ A (θ)} is that they correspond to the variational extremum (i.e. the saddle-point) of the free-energy functional for a specific choice of chemical potentials. Employing quasiparticle densities as variational objects and eliminating the variation of hole densities, using δρ A = −K AB δρ B (cf. Eq. (62)), the equilibrium partition integral is specified by the free energy functional Here each quasiparticle mode in the spectrum has been assigned a chemical potential µ A (θ). Specifically for the grand-canonical Gibbs equilibrium states, these have the form which can be easily inferred by resolving the energy density in terms quasi-particle modes, lim L→∞ (E/L) = e A ρ A , where e A (θ) denote bare quasiparticle energies carrying bare Cartan charges q (i) A . For the class of Hamiltonians (45) (with the anti-ferromagnetic exchange coupling), we find with In the large L limit, the partition integral (70) called the canonical TBA equations. As customary, we have introduced here a set of new thermodynamic variables, which we call the thermodynamic Y-functions; they play the role of 'Boltzmann weights', and thus in canonical Gibbs equilibrium assume the form where we have introduced the dressed quasiparticle dispersions with quantities ε A (θ) and q (i)dr A corresponding to the physical energy and the so-called dressed charges of the quasiparticle mode A, respectively. Accordingly, one defined the Fermi occupation functions n A (θ) as ratios Microscopic details of the model under consideration enter into equations (62) and (75) implicitly through the scattering kernels K AB (θ). We shall demonstrate in turn that the underlying Yangian symmetry Y (g) can be explicitly exhibited by inverting the integral kernel, namely by computing the Fredholm resolvent R AB (θ), defined through where the identity operator on the right-hand side is understood as 1 ≡ δ AB δ(θ).
Unitary spin chains. For concreteness, we consider here explicitly the unitary cases g = su(n). For the basis of Weyl generators ascribed to simple roots we can take E αa = |a a + 1|, such that ath Cartan element H αa takes the form With this choice, the Killing metric κ coincides with the Cartan matrix, κ ≡ K su(n) . The Fredholm resolvent has a remarkably simple structure where the s-kernel (defined as the solution to equation K 1 − s K 2 = s) reads explicitly and I AB ≡ I (a,s),(a s ) = δ a,a I A∞ s,s + δ s,s I is an adjacency tensor of the form A n−1 × A ∞ , where I An pertain to incidence matrices of dimension n − 1 associated to A n -type Dynkin diagrams whose 'bulk part' reads I s,s = δ s,s −1 + δ s,s +1 . The transformed source term are given by functions As explained in the subsequent section, asymptotic quantities y (∞) a depend solely on the U(1) chemical potentials h i assigned to the Cartan charges Q (αi) (associated with ith simple root α i ).
For example, in the grand-canonical Gibbs ensemble in a homogeneous SU(n) quantum chain with onsite representations V spωa p , the local source terms read explicitly d a,s (θ) δ a,ap δ s,sp β s(θ), modulo a multiplicative factor which depends on normalization.

Universal dressing equations
We would now like to present a more suggestive physical interpretation of the outlined TBA formalism, revealing itself upon interpreting Eqs. (62) and (75) as an nonperturbative renormalization of the quasiparticle's bare dispersion laws. Transforming quasiparticle's bare properties into the dressed ones effectively amounts to taking into account the net contribution from macroscopically many elastic two-body interactions with a thermal background. To be more specific, to every bare charge density q A (θ) per quasi-particle mode A = (a, s) (with rapidity θ) we associated the dressed counterpart q dr A (θ) by applying the 'dressing transformation' It is remarkable that such a dressing transformation is a particular functional of the equilibrium Fermi function which directly reflects the structure of the underlying symmetry algebra. The first thing to note is that the Bethe-Yang integral equations (62) for the total densities of available states and the TBA integral equations (75) for the equilibrium free energy are indeed merely two different manifestations of one and the same dressing transformation: the former one in fact encodes information about the dressed momenta via while the latter one provides the dressed quasi-particle energies. Indeed, the total state densities precisely coincide with dressed momentum derivatives, which is transparent already at the level of Eq. (53) (being the finite-volume counterpart of Eq. (91)). The dressing operator therefore admits the following compact universal form Dressing of Cartan charges. The dressed values of the Cartan charges can be most easily inferred from the thermodynamic Y-functions using Eq. (77), that is Several remarks are in order at this point. We first wish to stress that (i) the dressing operator (92) does not in general commute with rapidity differentiation and that (ii) it should not be confused with conventional dressing of excitations provided by the so-called shift function (or backflow), cf. [87,121]. Let us remind that the presence of finite chemical potentials h i > 0 breaks G-invariance of a state down to U(1) ×r generated by the Cartan elements. To uniquely specification elementary quasiparticles, one has to select a basis of simple roots; there are different choices which are related by a discrete symmetry transformation of the root system called Weyl group. For instance, in the case of su(n) the latter is just the permutation group Unitary quantum chains. In analogy to Eqs. (86), the dressing transformation for local charge densities q A (θ) can be also be cast in a concise group-theoretic form. To this end, we consider below the simplest case of simplylaced algebras su(n). Introducing the 'Baxterized' Cartan matrices, the dressing transformation can presented in the form where for the moment F a,s are just some dummy functions. One can easily verify that equations (62) correspond to the special case and can be retrieved by setting Employing for convenience rapidity derivatives, q A (θ) ≡ ∂ θ q A (θ) ≡, the dressing transformation can be brought into the following compact universal form [44] K dr where we have simultaneously introduced the following dressed (i.e. state-dependent) analogues of the Baxter-Cartan matrices K dr s,s [n a,s ] ≡ δ s,s − s(θ)I A∞ s,s n a,s , K dr a,a [n a,s ] ≡ δ a,a − s(θ)I An a,a n a ,s .
Spin-S SU(2) ferromagnetic chains. To demonstrate usefulness of the above formulation, we shortly consider as a concrete example the structure for the case of SU(2) quantum chains of spin-S degrees of freedom (with 2S ∈ N). Since r = 1, there is a single magnon species in the spectrum. Let us moreover denote by p (2S) s (θ) bare momenta of bound states (s-strings) in a spin-S ferromagnetic state, and furthermore introduce for convenience the 'G-tensor' obeying the matrix convolution identity (1 + K) G = s. The above F -functions can be identified with their dressed counterparts, that is F (2S) s ≡ G dr s,2S . By applying the dressing operator, one arrives at the system of coupled equations [60] s −1 δ s,s − I A∞ s,s n s F where the source term in Eqs. (101) 'jumped' at the 2S-th node owing to To solve Eqs. (101), the strategy is to first solve the homogeneous part of the recurrence relation and then solve the 'gluing condition' at the 2S-th node. Explicit results can be found in Ref. [60] and we do not reproduce them here. Let us just note that, thanks to the kernel identity K dr s,s = G dr s,s −1 + G dr s,s +1 , the complete set of solutions F (2S) s for all 2S ∈ N allow to reconstruct the complete information about the dressed quasiparticle scattering kernels K AB in the fundamental (Heisenberg S = 1/2) SU(2) chain [60].

E. Infinite temperature limit
Every admissible set of quasi-particle densities ρ A (θ) bijectively corresponds to thermodynamic Y-functions Y A (θ). State-specific information thus reveals itself upon analyzing the analytic structure of the Y-functions. To this end Y A (θ) are analytically continued them into the complex domain (i.e. θ-plane), where they exhibit isolated zeros and poles (or even branch cuts [120], in principle). This analytic data is induced by the structure of the local source terms source terms d a,s (θ), which are, in turn, directly linked to values of chemical potentials µ A (θ) specifying a (generalized) Gibbs ensemble.
It is rather unfortunate that even in the presence of integrability, the TBA equations (75) do not in general permit non-perturbative analytic solutions. At this junction one is typically resorts to numerical evaluations and solves a truncated system of coupled integral equations with use of an appropriate iterative scheme. There is nevertheless an important exception in this regard, corresponding to grand-canonical Gibbs states in the limit of infinite temperature β → 0. In lattice models invariant under a non-abelian compact Lie group G such a state is always well defined. Since as β → 0 the source terms evidently disappear, d A → 0, the thermodynamic Y-functions Y A (θ) become accordingly flat (i.e. constant functions of rapidity θ) and, to distinguished them explicitly from Y-functions, we denote them subsequently by Y A . In effect, the TBA dressing equations severely simplify and take the form of algebraic equations which can be solved in closed form. The remainder of this section is devoted to studying this particular case.
Unitary series. Let us briefly return to the basic example of the classical unitary series su(n). Using property 1 s = 1 2 , Eq. (86) can be brought in the form where for clarify we have suppressed dependence on the r = n − 1 U(1) chemical potentials. Algebraic relations (103) are the celebrated Y-system relation (originally the context of certain integrable QFTs [122], see also [123]), presently specialized for the absence of the the spectral parameter θ. While finite temperature and in general the thermodynamic Y-functions inherit non-trivial θ-dependence through the local source terms d A (θ), this nonetheless does not modify their large-s asymptotics, which is fully determined by the set of distinguished U(1) chemical potentials h i (for i ∈ I r ). This particular property makes the task of studying the large-s asymptotic scaling properties much easier, enabling to carry out the entire analysis in the algebraic setting of the infinite-temperature grand-canonical Gibbs state (given by Eq. (103)). For the particular case of su(n) algebras, the solution Y a,s ({x i }) of Eqs. (103) can be expressed as ratios of classical rectangular characters χ a,s (g 0 ), with χ 0,s = χ n,0 = 1, where diag(x 1 , x 2 , . . . , x n ) (with n i=1 x i = 1) parametrizing the maximal torus T , i.e. the maximal abelian subgroup of G = SU(n). Classical su(n) characters obey the constant Hirota bilinear relation [89,90,104] χ 2 a,s = χ a,s−1 χ a,s+1 + χ a−1,s χ a+1,s .
We provide additional details, including explicit determinant expressions, in appendix C 2. As a corollary, the infinite-temperature Fermi occupation functions admit a simple algebraic form

Yangian characters
In order to further examine the formal structure of the TBA Y -functions, we make a short excursion into the theory of representations of Yangians Y (g). As already outlined previously for the particular case of su(n), it is remarkable that the Y -functions Y A (and the associated Fermi functions n A ) associated to Yangian symmetry Y (g) for any simple Lie algebra g can be expressed in terms of classical g-characters [80].
Thermodynamic Y -functions Y a,s can be expressed as certain nonlinear transformations of T -functions as follows An infinite set of T -functions T a,s here formally correspond to characters associated with rectangular representations of Yangians Y (g), namely T a,s = χ(W a,s ), satisfying a variant of the Hirota bilinear relation [83] T 2 a,s = T a,s−1 T a,s+1 + T a,s .
For compactness of notation, we keep dependence on chemical potentials x i of the Cartan U(1) charges implicit. The last term in Eq. (110) can be formally expressed as for an appropriate g-dependent 'adjacency tensor' A which can be found in e.g. [83] (there, and elsewhere in the literature, Eq. (110) is referred to as the 'unrestricted Q-system', which however should not be confused with functional relations amongst the Baxter Q-functions).
In the simplest case of simply-laced g, the last term in Eq. (110) assume a simple form T a,s = b∼a T b,s , with ∼ running over all the adjacent nodes in the Dynkin diagram. The exact form of T -systems functional relations (110) are spelled out in appendix C 4. There we explain how T a,s can be expanded in terms of g-characters χ a,s for the case of classical Lie algebras g, and exemplify it on a number of simple cases.
Remark. It is instructive to mention that functional relations (110) are just a special case of a more general system of functional relations for rapiditydependent thermodynamic T -functions T a,s (θ), namely the θ-dependent counterpart of Eq. (110). The latter go commonly under the name of the T -system [124] (see [107] for an overview), which are omnipresent in the literature of exactly solvable models of statistical mechanics. Thermodynamic state functions T a,s (θ) can be understood as a distinguished gauge-covariant parametrization of the thermodynamic Y-functions. Every vacuum solution to the T -system relations represents a distinct equilibrium macroscate, where state-dependence is reflected through analytic properties of T -functions in the 'physical strip' P = {θ ∈ C; |Im(θ)| < 1 2 }, with large-θ asymptotics lim |θ|→∞ T a,s (θ) = T a,s . From a more physical perspective, T a,s (θ) can be identified with eigenvalues of (inhomogeneous) commuting column transfer matrices which associated to the 'auxiliary dimension' of a 2D vertex-model realization of a GGE density matrix (see Ref. [120] for details and explicit construction). It is thus important to not confuse functions T a,s (θ) with eigenvalues T a,s (θ) of commuting row transfer matrices T a,s (θ) which act on the physical Hilbert space.

F. Asymptotic analysis
After having place the TBA formalism firmly in place, we now return back to superdiffusive transport. We remind that the key properties determining charge transport are the asymptotic behavior of the TBA functions for large strings s → ∞. To this end, we subsequently examine certain formal properties of the dressing equations for the specific algebraic case of the grand-canonical Gibbs equilibrium ensembles in the limit of infinite temperature. The subject of study are states with unbroken symmetry G, which requires to set all the U(1) chemical potentials to zero, that is x i → 1.
To begin with, we notice that T -functions T a,s are, by construction, polynomials in s (for all a ∈ I r ). This readily implies that the Y -functions Y a,s are rational functions which only depend on quantum numbers (a, s). By virtue of Eqs. (107), the Fermi functions scale as n a,s ∼ 1/s 2 (and hence Y a,s ∼ s 2 ) at large s. This is a direct corollary by the Hirota relation (explicit formulae are given in appendix C 3). This in turn implies the asymptotic relation ρ tot a,s ∼ s 2 ρ a,s . While both state densities ρ a,s (θ) and ρ tot a,s (θ) depend algebraically on s, in distinction to Fermi functions they admit also non-trivial dependence on θ (even in the infinite temperature limit).
Equipped with this knowledge, we are now in a position to analyze the large-s behaviour of the dressing equations. We first inspect the basic case of the spin-1/2 Heisenberg chain SU(2) with a single type of magnon species. It will be sufficient to analyze for the canonical form of the Bethe-Yang equations Two observations can be readily made: (i) for fixed value of θ K s (θ) = 2π p s (θ) ∼ 1/s, and (ii) due to telescopic cancellation of poles the differential scattering phase shifts K s,s involve 2 min(s, s ) terms. Based on this, we can make the following estimate at large s ρ tot where each term under the sum has been upper-bounded by a constant using that 1 K s = 1. The second sum can be understood as the residual contribution coming from large s -strings with s > s, which gets suppressed in the ∼ 1/s fashion. The first sum converges in the s → ∞ limit provided lim s→0 ρ tot s = 0, i.e. the total density of states ρ tot s has to decay to zero at large s in an algebraic fashion. Using finally that the bulk recurrence relations pertaining to the algebraic dressing equations involve only rational functions of s, this implies ρ tot s ∼ 1/s scaling at large s.
In more general integrable models (with nested spectra), one has to additionally account for the quasiparticles of different flavors, whose mutual interaction are prescribed by the Cartan matrix K. The above argument can be easily adapted by including an extra sum over the flavor index.

G. Analytic solutions
To solidify the conclusion of the preceding section, we now work out the exact closed-form solution to the algebraic equations. We note that the basic case of the SU(2) spin-S (Babujian-Takhtajan) quantum chains, the full solution has already been obtained in ref. [60]. Our primary interest here are model that possess Lie symmetries of higher rank.
For definiteness, we specialize below to SU(n) quantum chains with fundamental onsite degrees of freedom. While other cases can be treated in a similar fashion, for compactness of presentation we postpone a complete and comprehensive analysis for a separate technical study.
Let us consider, as our starting point, the algebraic limit of the group-theoretic dressing equations (97). In Fourier k-space, these turn into a system of coupled algebraic equations of the form with s −1 (k) = 2 cosh (k/2). The first step is to solve the homogeneous part, which can be achieved with the following ansatz F a,s (k) = A a,s (k)K a+s−1 (k)+B a,s (k)K a+s+1 (k), (115) where K s (k) = e −s|k|/2 are the elementary scattering kernels in Fourier space. Plugging the ansatz into Eqs. (114), Laplace transforming from the k-plane to the z-plane, and demanding the null condition for all the residues located at z i ∈ 1 2 N, we end up with following system of recurrence relations A a,s =n a,s−1 A a,s−1 + n a−1,s A a−1,s , (116) B a,s =n a,s+1 B a,s+1 + n a+1,s B a+1,s , (117) A a,s − B a,s =n a,s+1 A a,s−1 + n a+1,s A a+1,s −n a,s−1 B a,s−1 + n a−1,s B a−1,s , where the Fermi functions are given by su(n) characters as per Eqs. (107). The solution has the form A a,s = A 0 χ a,s χ a−1,s−1 χ a−1,s χ a,s−1 , B a,s = B 0 χ a,s χ a+1,s+1 χ a+1,s χ a,s+1 , for some yet to determined coefficients A 0 and B 0 . By plugging these back to Eqs. (118), we retrieve the Hirota relation for classical characters (106). To fix the undetermined constants, we find the particular solution of Eqs. (114), with the source term attached at the first node at (a, s) = (1, 1), yielding Transforming back to the θ-plane, we finally find F a,s (θ) = χ a,s χ 1,1 χ a−1,s−1 χ a−1,s χ a,s−1 K a+s−1 (θ) − χ a+1,s+1 χ a+1,s χ a,s+1 K a+s+1 (θ) .
At the end we are interested in the x i → 1 limit, in which the su(n) characters become dimensions d a,s or rectangular irreducible representations V s ωa (see appendix C 2 for details). Specifically, the total state densities read ρ tot a,s (θ) = s + n − a n There is an analogous expression (up to a multiplicative prefactor) for the dressed differential quasiparticles energies ε a,s (θ) which can be obtained by replacing K s (θ) with K s (θ). In rapidity space, the latter is given by a discrete family of Cauchy distributions Notice that any finite fixed rapidity θ, the density of states ρ tot a,s decays as ∼ 1/s 2 in the large-s limit. However, when integrated against any dummy integrable function f (θ) as in GHD, we have the following large-s properties for every quasiparticle flavor a.
As announced in Sec. II C, these properties allow us to infer the superdiffusive scaling of charge dynamics with an algebraic dynamical exponent z = 3/2.
Non-fundamental onsite representations. The large-s scaling properties (125) and (126) likewise hold in models with non-fundamental local degrees of freedom, that is for any on-site, finite-dimensional unitary irreducible representations V Λ . At the level of algebraic dressing equations, this amounts to moving the source terms to a generic position. For a quick illustration, we derive the explicit solution to the SU(3) chain in the adjoint representation (8), with Dynkin labels Λ p = ω 1 + ω 2 ≡ [1, 1]. By virtue of Z 2 invariance of the dressing equations (under interchanging the flavors) we have F 1,s = F 2,s , and therefore where I A∞ s,s = δ s,s −1 + δ s,s +1 . The solution reads V. SEMICLASSICAL SPECTRUM Thus far we have established that the anomalous charge transport comes from giant quasiparticles (i.e. bound states carrying large quanta s) immersed in a charge-neutral thermal background. To give a different angle and better elucidate their physical nature, we next discuss a semiclassical description of these giant quasiparticles. It is helpful to first review some results along these lines that were shown for SO(3) spin chains [42]. We then proceed by outlining how these arguments generalize to other symmetry groups.
We begin with a classical SO(3) ferromagnet in its ground state, i.e., all spins aligned in some direction. The elementary excitations of the ferromagnet are Goldstone modes in which the spin orientation changes slowly across the lattice. One can regard any sufficiently slowly varying spin texture as being composed of Goldstone modes, as it locally consists of slow modulations (rotations) of the vacuum orientation. In one dimension, the longwavelength dynamics of these ferromagnetic, quadratically dispersing, Goldstone modes is governed by the Landau-Lifshitz equation, which is nonlinear but integrable PDE possessing stable soliton solutions of any width w [125]. The properties of these soliton solutions can be inferred by observing that they are wavepackets of Goldstone modes that are stabilized by nonlinearity. For example, a soliton of width w is made up of Goldstone modes with characteristic momentum 1/w and thus characteristic energy density 1/w 2 ; therefore its characteristic energy scales as 1/w and its velocity also as 1/w. On the other hand, its spin orientation is away by O(1) from the vacuum so it carries spin ∼ w. These "bare" properties of classical solitons precisely match those of the large-s Bethe ansatz strings, suggesting that there is some correspondence between them. Indeed, the exact correspondence between large-s strings and large-w solitons can be explicitly established, not only at the level of individual spin-wave configurations as usual, but even in high-temperature equilibrium ensembles [42].
This classical quasiparticle construction generalizes to arbitrary quantum or classical ferromagnets. As discussed below, a distinguished feature of such ferromagnets are quadratically dispersing Goldstone modes. One can analogously construct slowly varying wavepackets of these Goldstone modes; by construction these have the same scaling as the large-s strings, so it is tempting to identify the two types of excitations in the general case also. Although we have not yet attempted to analytically solve the classical problem in full generality, by noting that large-s strings must correspond to objects in the semiclassical spectrum and that Goldstone-like excitations exhaust this spectrum, it is tempting to infer the existence of such solitons from the integrability of the quantum model. However, one major puzzle arises when one tries to perform such an identification: the number of magnon flavors (which is always set by the rank r of the group G) is far fewer than that of Goldstone modes; the latter is given by one half of the real dimension of a coset manifolds G/H, where H is the residual symmetry group of the ferromagnetic state. The "missing" Goldstone bosons must in some fashion emerge out of the Bethe ansatz spectrum; however, they cannot be identified directly with elemen-tary (physical and auxiliary) magnons of integrable quantum chains. Curiously, they rather correspond to 'stacks' of magnons, which are long-lived wavepackets in a finite system, and only become sharp eigenstates in the infinite system limit.
The main result of this section is identifying the correspondence between Goldstone modes in the semiclassical spectrum and stacks of magnons (and condensates thereof, representing semi-classical Bethe s-strings) in the Bethe ansatz spectrum. Using this correspondence, one can identify large-spin semiclassical excitations above the ferromagnetic vacuum precisely as stacked strings. These semiclassical excitations can then be identified as solitons that involve smooth maps from R to the coset space G/H. While it would be interesting to construct such solitons directly in the classical theory, this task is left to future work.

A. Goldstone modes
As a natural starting point to describe classical continuous ferromagnets, we being by characterizing the spectrum of linear fluctuations of a ferromagnetic order parameter. This leads to the notion of Goldstone excitations which naturally arise in systems with spontaneously broken continuous internal symmetry. To this end, we begin by shortly revisiting the Goldstone theorem in a general context of non-relativistic field theories. Our objective here is mainly to characterize the (non-relativistic) Goldstone modes that govern the low-energy spectrum for the class of ferromagnetic quantum spin chains introduced in Sec. IV B.
First we recapitulate the main result on the counting of Goldstone modes systems without Lorentz invariance, independently worked out in Refs. [126,127] and [128] (optimizing the earlier Nielsen-Chadha inequality [129]). The conventional setting concerns Hamiltonian systems invariant under a non-abelian Lie group G, possessing a degenerate ferromagnetic ground state manifold. Spontaneous breaking of internal symmetry amounts to picking a particular vacuum polarization (order parameter), say Ω, thereby breaking the symmetry of the isometry group G down to stability subgroup H ⊂ G of Ω, determined by condition We are specifically interested in homogeneous (i.e. translational invariant) ferromagnetic chains (of length L), where the global (pseudo)vacuum is simply a product state Ω L = L =1 Ω. It is therefore sufficient to carry out the analysis at the level of local Hilbert spaces and (with no loss of generality) we can assume Ω ≡ |0 0|.
At the algebraic level, the symmetry breaking pattern G → H implies the splitting g = h ⊕ m, with h being the Lie algebra that generates H (which contains the Cartan subalgebra t) and and m denoting the linear span of the 'broken generators' (these do not enclose an algebra). The number of broken generators is Let X σ ∈ m, σ ∈ {1, 2, . . . , dim(G/H)} be a hermitian basis of m. By the Goldstone counting theorem, the number of independent Goldstone modes n GB is with matrix storing the vacuum expectation values of the commutators amongst the broken generators X σ . Goldstone modes come in two varieties, depending on the type of the dispersion law at long wave-lengths (k → 0), ω(k) ∼ k ν : type-I (or type-A) with ν odd, and type-II (or type-B) with ν even. In general, we thus have with In relativistic systems, Lorentz invariance forces V to be identically zero. In contrast, target spaces of nonrelativistic ferromagnets are symplectic manifolds. These have even (real) dimension and consequently V has full rank. Therefore n I = 0 and In simple terms, every canonically-conjugate pair of broken generators contribute an independent magnon branch, a quadratically dispersing long-wavelength spinwave excitation of the ferromagnetic order parameter.

B. Symmetry breaking patterns
We now describe the spectra of Goldstone modes for the class of G-invariant ferromagnetic chains. The structure of coset spaces depends on the choice of a finitedimensional irreducible onsite representation V Λ . Let us for simplicity suppose for the moment that the local degrees of freedom transform in the defining representation of g. The highest weight is the first fundamental weight Λ = ω 1 = [1, 0, . . . , 0] of g, and the stabilizer H ⊂ G of the vacuum state Ω has the structure U(1) × H (i.e. H is not semi-simple). For the classical series, the coset structure is summarized in Table I. Next, we specialize to the family of fundamental representations with Λ = ω a . In this case, the symmetry breaking pattern can be inferred directly from the enumerated Dynkin diagrams with use of a graphical recipe by simply breaking the bonds which connect to the node containing the non-zero label. This is illustrated in Fig. 5 on a two simple examples. The list of possibilities for all the classical series is summarized in Table II. Employing the same trick allows us to determine the structure of coset spaces also for G-invariant classical ferromagnets corresponding to the exceptional groups G 2 , F 4 , E 6−8 , specializing to the fundamental representations of g (for completeness we include them in appendix D). The outlined recipe likewise applies for generic (nonfundamental) finite dimensional irreducible representations V Λ of g with Dynkin labels Λ = [ω 1 , . . . , ω r ]; the rule is to break all the bonds that connect to nodes with nonzero Dynkin labels. The resulting coset spaces are generalized flag manifolds F i1,i2,...,i k ≡ G/H i1,i2,...i k , where the indices i i mark the nodes of the Dynkin diagram of g with non-vanishing Dynkin labels (that is i i = 1 if and only if m i = 0). Compact complex manifolds F i1,i2,...,i k indeed exhaust all Kähler manifolds (i.e. possess compatible Riemannian metric and symplectic structures) which are homogeneous spaces with a transitive action of G. Furthermore, points on such flag manifolds are in oneto-one correspondence with generalized coherent states in a specific representation of group G, namely for the highest weight state |Λ of V Λ . If some m a = 0 however, not all coordinates w − α are independent and some can be eliminated. For instance, in the unitary case (A-series), the family of generalized flag manifolds has the structure SU (n)/S(U (n 1 ) × U (n 2 ) × · · · × U (n j ) × U (1) ×k ), (137) for integers j, k ≥ 1 and n 1 ≥ n 2 ≥ . . . ≥ n j > 1, subjected to k +  (6)) or H 1,1 = U (1) ×2 (as e.g. for (8), (15) or (27)). The upshot of that is that physical degrees of freedom of the lowenergy effective theories corresponding to quantum spin chains whose degrees of freedom transform in different irreducible representations can live on the same target coset manifold. Below we explain this fact from another perspective by directly examining the structure of semiclassical spectrum.

C. Counting degrees of freedom
Magnon modes can arise as plane-wave solutions to linear partial differential equations. On the other hand, classical field theories are interacting systems governed by nonlinear equations of motion which possess a vastly richer spectrum of solutions. A hallmark feature of nonlinear completely integrable PDEs are soliton solutions. Solitons refer to nonlinear ballistically propagating field configurations that behave as particles, whose stability is ensured by integrability. In isotropic ferromagnets considered in this work, magnetic solitons assume a nontrivial internal structure. Indeed, from the perspective of quantum spin chains, soliton emerge as certain macroscopic coherent superpositions of magnons. Such states represent highly-excited semiclassical eigenstates in the the low-energy spectrum of the model and their time evolution is generated and an effective classical Hamiltonian. This fact can be established in various ways; the most standard and direct is employ a standard pathintegral technique to project the quantum many-body Hamiltonian onto the manifold of (generalized) coherent states |ψ , given formally by the 'vacuum rotation' |ψ = g |Λ , as prescribed by Eq. (136). This way one can deduce a classical evolution for the ferromagnetic order parameter that takes values on a quotient manifolds G/H, where H is the vacuum stability subgroup defined above in Eq. (129) (see also e.g. [40,130,131]). Classical fields can be thus identified with coordinates of generalized (Perelomov) coherent states [132].
Every admissible classical solution is a particular longwavelength (i.e. low-momentum) solution of the Bethe Ansatz equations with a macroscopically large number of quanta. This particular regime, often referred to as the 'asymptotic Bethe Ansatz', has been first studied by Sutherland [46] and subsequently more thoroughly examined in [133]. Nevertheless, an exact identification with classical field solutions has however only been made precise afterwards in Ref. [134] which provided a prescription to construct the associated classical spectral curve (the characteristic equation associated with the classical monodromy matrix, see e.g. [135,136]) by transcribing the asymptotic Bethe equations as a Riemann-Hilbert problem. This method permits to describe and classify the nonlinear modes for the class of so-called finite-gap solutions [137].
In the remainder of this section we explain how the spectrum of linear and nonlinear mode of classical integrable ferromagnets arises from the nested magnon spectra in integrable G-invariant ferromagnetic quantum chains. Complete characterization of classical algebraic curves that emerge in this picture is a technical matter exceeds our scope and we thus do not undertaking it here. We rather focus here on a deceptively innocent problem of identifying the relevant emergent classical degrees of freedom. The task turns out to be quite subtle and in order to resolve we will need a basic description of the asymptotic Bethe Ansatz.

Asymptotic Bethe Ansatz
We begin by recalling that every finite-volume (highest-weight) Bethe eigenstates in a quantum chain can be uniquely resolved in terms of magnons excitations. Individual unbond magnons correspond to realvalued Bethe roots with rapidities θ ∼ O(L 0 ). As described earlier, complex-valued rapidities of constituent magnons which are a part of a bound state take the form of Bethe strings, see Eq. (59). Semiclassical eigenstates in the spectrum belong to highly-excited states with O(L) excitations and large rapidities θ ∼ O(L), which can be thought of condensing into a finite number of coherent macroscopic modes. In this low-energy regime, magnon momenta are inversely proportional to the system size, p(θ) ∼ O(1/L). By accordingly introducing rescaled rapidities λ As usual, n (a) j represent (integer) mode numbers due to multivaluedness of complex logarithm.
Classical nonlinear modes pertain to bound states of N a ∼ O(L) constituent low-momentum magons that all share the same mode number. Crucially, rapidities of nearby magnons are separated by an amount o(L); however, since θ ∼ O(L) magnons need not be spaced equidistantly as the Bethe strings with θ ∼ O(1). The situation one encounters is in general as follows: upon taking the large-L limit, rapidities λ (a) j distribute densely along certain one-dimensional disjoint segments (i.e. contours) C (a) = i C (a) i in the spectral λ-plane, described by some smooth (in general non-uniform) densities a (λ) of Bethe roots λ (a) j supported on C (a) . Such macroscopic low-energy configurations assume a purely classical description as they correspond to interacting 'nonlinear Fourier modes' of nonlinear classical field configurations. The long-wavelength spectrum of excitations involves in addition non-macroscopic excitations, i.e. modes that carry an infinitesimal amount of conserved (classical) charges -these are the aforementioned Goldstone modes (i.e. fluctuations above the ferromagnetic vacuum) and correspond to infinitesimally small contours (isolated poles) in the spectral complex plane.

Classical degrees of freedom
Based on the above picture, it is tempting to readily conclude that the total number of internal (isotropic) degrees of freedom of a generic classical field configurations coincides with the number of distinct flavors of the Bethe roots (which equals r = rank(g)). Before long one however realizes why such a naïve identification is flawed. Here we mention two apparent inconsistencies that arise: 1. Imagine first, for definiteness, the fundamental SU(n) ferromagnets. In this case the classical low-energy theories are given by the higher-rank Landau-Lifshitz equation constrained to complex projective manifolds CP n−1 ∼ = SU(n)/[U(1) × SU(n − 1)] (the derivation can be found in e.g. [40]), with classical phase space of real dimensions dim(CP n−1 ) = 2(n − 1). The Goldstone theorem predicts n − 1 magnonic branches in this case, and indeed this number precisely agrees with rank of SU(n) which also matches the total number of flavors in the SU(n) magnets. Although this appears to be a good sign, we nevertheless run into a problem: in generic excited quantum eigenstates, auxiliary quasiparticles cannot be excited on their own without first exiting physical (i.e. momentum-carrying) excitations. This can also be readily deduced from the selection rule (58).
2. Consider next, as another example, the SO(5) chain of rank r = 2 in the fundamental onsite representation V . The low-energy is associated with the coset space SO(5)/[SO(2) × SO (3)] of real dimension 6. This time we expect to find 3 Goldstone modes which this time obviously exceeds the total number of flavors. More generally, one would expect the number of distinct types of macroscopic semiclassical strings to always equal r = rank(g), which is typically vastly smaller than the number of independent Goldstone modes (cf. Sec. V B).
The only viable option to avoid an apparent paradox is to look for additional classical degrees of freedom that somehow managed to stay unnoticed. These clearly would have to be emergent as our enumeration of the genuinely quantum excitations spectrum is unquestionably correct. This brings us to the notion of stacks. The existence of stack excitations has, to the best our knowledge, been first encountered in Refs. [138][139][140] in the context of gauge-string duality, devoted specifically to the lowenergy sector for the particular case of unitary supergroups related to the underlying symmetry algebra of the superconformal Yang-Mills theory [141,142]. Outside of that, it appears there a general or systematic discussion devoted to these emergent inherently classical objects is still lacking in the literature. To this end, we proceed to explain how the notion of stacks naturally arises in the context of integrable quantum ferromagnets with continuously degenerate ground states.
In most general terms, stacks pertain to excitations with an internal multi-flavor structure. They are produced by merging together a subset of Bethe roots of distinct flavors that share the same rapidity (no merely just the real part). As already emphasized, in a quantum spin chain such a configuration is not permitted as rapidities must be pair-wise distinct. In the low-momentum scaling limit described by the Asymptotic Bethe Ansatz this restriction is however lifted as Bethe roots of different flavors may approach arbitrary close to each other before eventually recombining into an independent lowmomentum fluctuation mode. Although a large number of stacks can form this way, only those stacks that involve momentum-carrying roots can be regarded as physical fluctuations, whereas stacks made of solely from auxiliary excitations are likewise called auxiliary. Similarly to the elementary excitations of a quantum chain, stacks can exert attractive interactions and therefore can also condense into 'large' (i.e. non-linear) modes which carry finite filling fractions and charge densities. For the complete classification of classical modes additional quantum numbers are thus required to properly account for the internal structure all the admissible stacks. This additional information, which can be though of various polarization directions in which a classical configuration can vibrate, depends on both the symmetry algebra and on-site representation under consideration. For illustration, let us have a look at our main example of the fundamental SU(n) chains. The asymptotic Bethe Ansatz equations, see Eq. (138), take place on an n-sheeted Riemann surface, whose sheets are stitched together at a finite number of branch cuts (see e.g. [134,135] for ad-ditional details). Standard branch cuts, forming along contours C (a) i along which the Bethe roots condense, are of the square-root type. Creating a single excitation of flavor a amount to connects the ath and (a + 1)th sheet with an isolated pole (infinitesimal cut). Exciting a stack (a, b) means joining two non-adjacent sheets indexed by a and b + 1. Equations (138) need to be appropriately amended to incorporate these extra stack excitations. In Fig. 6 we illustrate the process of stack formation on a four-sheeted Riemann surface associated to the asymptotic Bethe equations of the SU(4) quantum chain made out of fundamental spins (4).
In summary, the number of independent classical degrees of freedom equals the number of distinct physical stacks. In the next section, we carry out a systematic analysis to determine the admissible symmetry breaking patterns and classify the corresponding physical stacks.
Remark. We would like to re-emphasize that stack excitations are are not permitted in the spectrum of a quantum chain due to kinematic constraints. In other words, upon including quantum fluctuations (even perturbatively at the level of asymptotic equations (138)) stacks are immediately forced disintegrate into individual magnons.

Giant quasiparticles as classical solitons
We have thus far have not yet encountered or mentioned solition modes. As a matter fact, classical soliton modes are (strictly speaking) not a part of the nonlinear (finite-gap) spectrum. Nevertheless, they emerge as from elliptic waves in the limit of unit elliptic modulus, which invariably requires decompactification of the spatial dimension [134,143]. This process, when observed in the spectral λ-plane, corresponds to a singular event in which two adjacent density contours merge together, resulting in a perfectly vertically oriented contour with constant unit density of Bethe roots. Such object have attracted considerable attention and have most commonly been referred to as condensates, see e.g. Refs. [134,135,139,143]  (note however that these are distinct from ordinary condensates that form square-root branch cuts of finitegenus Riemann surfaces). From the algebraic curve viewpoint, such uniform condensates arise from coalescing two square-root branch cuts into a logarithmic branch cut. The giant quasiparticles from the preceding discussion, i.e. low-momentum bound state excitations with macroscopic number of quanta, are semiclassically quantized soliton modes.

Counting stacks
Now we return back to the 'counting problem'. What is left is to extract the correct number of Goldstone modes from the semiclassical spectrum of quantum ferromagnetic chains whose elementary (both physical and auxiliary) excitations carry r distinct flavor numbers (one per each simple root α i of g). One can generate different types of stacks simply by forming linear combinations of α i with non-negative integral coefficients. In other words, excitation operators associated with stacks correspond to positive roots α ∈ ∆ + (g). Among those, the physical ones are those that involve at least one simple root, associated with physical (i.e. momentum-carrying) excitations. In the fundamental models for example, that is only α 1 .
There indeed exist a nifty graphical representation for all the physical stacks in terms of partially ordered sets known as Hasse diagrams. Let us first outline the construction by specializing to fundamental irreducible onsite representations Λ = ω a . One starts by picking one of the bottom nodes (associated to simple roots) and proceeds by 'climbing' in the direction of the top node (being the highest root) by successively adding one simple root at a time. With this rule, several inequivalent paths through the Hasse diagram can be generated. Since each vertex is an independent Goldstone mode, their total number n G equals the total number of vertices visited along the way. More generally, that is in the case of non-fundamental irreducible onsite representations V Λ with Λ = r a=1 m a ω a involving multiple non-vanishing Dynkin labels m a , there will be multiple momentum-carrying magnon species in the spectrum. Now one has to simply superimpose all the sublattices that emanate from base nodes α i for all a ∈ I r for which m a = 0. The 'maximal' breaking occurs when none of the m a 's vanishes, as then the isotropy group H coincides with the maximal torus subgroup T of G. This scenario results in the maximal number of Goldstone modes n G = 1 2 dim(G/T ) = |∆ + |, rendering all stacks physical excitations.
Finally, in cases when g contains a node of connectivity 3 (that is types D r and E-family), the Hasse lattices extend into the third dimension (this is illustrate in Fig. 8 on the example of so(10)).

VI. INTEGRABLE QUANTUM FIELD THEORIES
Integrable ferromagnetic quantum chains discussed thus far do not exhaust the list of integrable quantum models which possess global symmetries of non-abelian Lie group. There are several other non-relativistic systems, such as most prominently integrable fermionic models that exhibit Lie supersymmetries (Z 2 -graded Lie algebras). It should therefore not come as a surprise that the phenomenon of charge superdiffusion with exponent z = 3/2 likewise occurs in there too (see [36,44] for concrete models). What is more remarkable however is that that anomalous charge dynamics also emerges in relativistic invariant integrable quantum field theories in two space-time dimensions, provided that elementary particles carry internal isotropic degrees of freedom (which   (4) and B3 = so(7) of rank 3: in α-basis, any positive root α ∈ ∆+ is decomposes as a integral linear combination of simple roots αi (level equals the sum of coefficients), whereas ω-basis provides an expansion in terms of the fundamental weights (with coefficients being the Dynkin labels).
classically takes values in G or coset manifolds G/H).
Even though Goldstone modes of Lorentz invariant systems always comprise of linearly dispersing (i.e. type-I) modes, one should keep in mind that we are interested here exclusively in the time-evolution of the Noether currents and not correlation functions amongst individual components of quantum fields. A crucial observation in this respect is that the temporal components of the Noether two-currents may be understood as (quantized) spin waves. Moreover, in finite temperature ensembles excitations of these internal interacting degrees of freedom undergo non-trivial dressing. The situation in fact mirrors that of the Heisenberg spin chains or its higherrank analogues we investigated earlier. We shall not attempt to give a comprehensive exposition but rather demonstrate the basic principles on the O(3) NLSM, representing a paradigmatic integrable interacting QFT with non-diagonal scattering. We shall briefly outline the differences in a few other IQFTs that display SU(2) symmetry, while postponing a detailed quantitative study of other integrable nonlinear sigma models on Riemannian symmetric spaces G/H [150] and their classical limits for another study.  (10)). There are 20 roots in total (black dots), arranged into 7 levels. The lattice extends into three dimensions as simple root α3 has connectivity 3. To find the α-pattern corresponding to a simple roots αi, one begins at αi and collects also the nodes by always proceeding forward in the direction of the highest root. The Z2 automorphism of the Dynkin graph reveals itself as the reflection symmetry with respect to the main diagonal of the cube passing though nodes 4 and 5.

O(3) nonlinear sigma model
To begin, we spell out some the main properties of the O(3) NLSM quantum field theory, employing the Hamitonian formalism. The quantum field n = (n x , n y , n z ), subjected to the nonlinear constraint n · n = 1, transforms in the vector representation of O(3). Since the vacuum stability subgroup w.r.t. polarization axis (say Introducing the momentum π conjugate to n and the angular momentum field m = n × π perpendicular to n ( m · n = 0), the Hamiltonian density (here without the topological Θ-term) has the form The conserved Noether two-current associated to global O(3) rotations is given by for µ ∈ {x, t}. The integrated magnetization density provides a local conservation law, (d/dt) dx m(x) = 0. The O(3) NLSM is well-known to arise as the effective low-energy theory of (Haldane) spin-S antiferromagnets, where one identifies v = 2J S and coupling strength is given by g = 2/S, assuming large S; the S 2 -fields n governs fluctuations of the staggered order parameter above the antiferromagnetic (Néel) ground state, whereas m pertains to ferromagnetic (spin wave) fluctuations. Elementary excitations of Eq. (139) are an SU(2) triplet of massive spinfull bosons with relativistic (bare) dispersion e(p) = p 2 + m 2 . As customary, relativistic dispersion of massive particles can be parametrized by a rapidity variable θ e(θ) = m cosh (θ), p(θ) = m sinh (θ).
Mass m gets generated dynamically through dimensional transmutation and vanishes in the large-S limit. At weak coupling (g → 0), the theory becomes effectively free with dim(G/H) = 2 massless (Goldstone) bosons.
The integrability of the model manifested itself though (i) infinitely many "hidden" conserved currents and (ii) a completely factorizable many-body S-matrix. Notice however that exchange of isotropic degrees of freedom upon elastic interparticle interactions renders the Smatrix non-diagonal, namely particles do not only pick up a U(1) phase but change their state (internal quantum numbers). This feature bears direct analogy with the nested spin chains discussed previously in Sec. IV C.
Algebraic diagonalization using Bethe Ansatz has been first carried out in Ref. [56] (see also [150,151] for a more comprehensive discussion). The usual trick is to resolve the exchange of isotropic degrees of freedom at expense of introducing additional, called auxiliary, spin wave excitations (magnons); the extended basic of excitations which renders the scattering diagonal. These auxiliary magnons are, despite being massless, regarded as independent quasiparticles. Imposing finite-volume boundary condition, one arrive at the following the NBA equations [56] belonging to the sector with N θ momentum-carrying physical excitations, N λ magnons (with associated rapidity λ). There is a single rational scattering amplitude reading S(θ) = (θ − iπ/2)/(θ + iπ/2). Equations (143) already make it transparent that interactions among magnons are identical to those in the isotropic Heisenberg SU(2) chain, the essential difference being that here spin excitations are pinned onto an underlying field (a kind of inhomogeneous background). Consequently, inter-particle interaction lead to formation of bound states (Bethe strings). Taking the standard thermodynamic limit (that is L → ∞ with ratios N θ /L and N λ /L finite) and introducing physical density ρ 0 (θ) and densities of auxiliary magnonic Bethe s-strings ρ s≥1 (θ), one finds following (decoupled) Bethe-Yang equations ρ tot s = δ s,2 s ρ 0 + s I A∞ s,s ρ s .
with the incidence matrix of D ∞ Dynkin diagram (see Fig. 9). Notice that the U(1) chemical potential h enters only implicitly via asymptotics of the Y-functions, lim s→∞ log 1 s log Y s = h. The O(3) model is fact belongs to an infinite family of integrable SU(n)/SO(n) sigma models with n − 1 massive nodes whose thermodynamic particle content gets arranged according to A n−1 × D ∞ Dynkin diagrams [150].
The limit of infinite temperature is no longer physically meaningful in IQFTs. Although there is no issue with computing the thermodynamic Y-functions in the β → 0 limit at half filling h = 0 (which are given by Y 0 = Y 1 = 2 and Y s≥2 = s 2 − 1), the issue arises when solving for the densities Eqs. (145) which yield divergent rapidity integrals in the UV regime. We purposely avoid imposing a momentum cut-off (e.g. as in Ref. [152]) as in general, it would break integrability and thereby spoil the late-time decay of charge correlations. To ensure convergence of rapidity integrals temperature β > 0 has to be taken into account in a non-perturbative way. Another possibility to overcome the issue is to consider instead integrable lattice regularizations, i.e. realizing IQFTs are certain continuum scaling limit of inhomogeneous (staggered) quantum spin chains (as exemplified for the case of SU(2) PCF in [147]).
In the scope of physics applications, it is more meaningful to perform a low-temperature expansion of Eqs. (145) and Eqs. (146). We refer the reader to [60] and only provide a few remarks here. In order to properly account for the effects of thermal fluctuations in a finite-density state (that means even at arbitrarily small temperatures T > 0), it is crucial to retain all the contributions of the spin-wave excitations. Magnonic excitations can be safely discarded only in the regime h/T 1 which can be well approximated by the semiclassical description developed in refs. [153][154][155]. In a close analogy to the Heisenberg spin chain, magnons propagating in an unmagnetized (or half-filled) finite-density background carry vanishingly small amount of dressed magnetization, and one accordingly finds a vanishing spin Drude weight with a diverging spin diffusion constant D s ∼ 1/h.
Topological term. It is well known that the O(3) NLSM admits an integrable deformation with the inclusion of the topological Θ-term with Θ = π [56,156] (describing the low energy limit of SU (2)-invariant spin chains with odd spin S; in general Θ = 2πS [157]). This appreciably modifies the spectrum as instead of a massive triplet one now finds an SU(2) doublet of massless elementary excitations with bare dispersion where here ± designate the right (p > 0) and left (p < 0) moving components, and M is an arbitrary mass scale (which is inessential as far as only the left-left and rightright scattering matters). At the level of TBA description, the physical species comprise of the left and right movers with densities ρ ± (θ) [56]. The Internal magnon structure however remains intact [56]. For instance, the TBA equations in the quasi-local form are now of the form For further details we again refer to ref. [60] and references therein. Classical picture. The origin of magnon excitations in relativistic systems can be exhibited already at the classical level. To this end we briefly consider the O(3) NLSM as an integrable classical field theory [158]. The Euler-Lagrange equation for the classical field n ∈ S 2 reads (in dimensionless units) [159] n tt − n xx + ( n 2 t − n 2 x ) n = 0.
Owing to Lorentz invariance, linear fluctuations above a degenerate ground state comprise of two transversal linearly-dispersing (type-I) Goldstone modes. We are however interested the time-evolution of the temporal component of the conserved Noether two-current which presently corresponds to the ferromagnetic order parameter m = n × π. Its equation of motion reads simply m t = n × π t . In terms of the Hamiltonian equations, one readily finds that magnetization field m satisfies thereby revealing the mechanism for how spatial fluctuations of n(x, t) generate the dynamics of the angular momentum field m(x, t). In the quantum O(3) NLSM, fluctuations of these magnetization waves carrying integer quanta scatter completely elastically of one another, as described by Eqs. (143).

Other integrable QFTs
The above basic example of the O(3) NLSM already displays all the characteristic features of IQFTs possessing isotropic degrees of freedom. This means that analogous types of TBA equations are found in other integrable QFTs with non-abelian isometry groups [149,160]. The best studied examples are O(n) NLSMs on O(n)/O(n − 1) ∼ = S n−1 target spaces, described by Lagrangians In the simply-laced cases O(2r) with r ≥ 4 [161], the thermodynamic spectrum of excitations comprises r flavors of auxiliary quasi-particles (one per node in the D r Dynkin diagram), each forming an infinite tower of magnonic bound states (s-strings). The O(4) NSLM (in the vector representation) is special as it can be viewed as SU(2) L × SU(2) R PCF with particles transforming in the bifundamental representation of su(2) [162]. Its spectrum involves massive spinfull particles with two types of SU(2) spins -these are elementary excitations above the Fermi sea (anti-ferromagnetic ground state) in an SU (2) spin-S chain. Its thermodynamic particle content is depicted in Fig. 9. From this perspective, the SU(2) PCF can be, loosely speaking, perceived as the QFT bosonic counterpart of the Fermi-Hubbard model [92,93]. In the strong-coupling limit, the SU(2) PCF model splits up into two copies of the isotropic Heisenberg chain. The classical limit and the associated Riemann-Hilbert equations can once again be derived from the asymptotic Bethe equations [160], governing the regime with N → ∞ particles with large rapidities θ ∼ O(N ) (with quantity m L = exp (−2π N ) playing the role of a small parameter). Another prominent class of IQFTs are SU(n) chiral Gross-Neveu models (cGN) [163]. Let us have a quick look at the simplest SU(2) case, describing two interacting Dirac fermions expressed in terms of two-component spinors ψ a (a = 1, 2) with Lagrangian density with / ∂ ≡ γ µ ∂ µ and γ-matrices γ 0 = σ 1 , γ 1 = iσ 2 , γ 5 = γ 0 γ 1 obeying thhe Clifford algebra γ µ , γ ν = 2η µν with metric η = diag(1, −1). Lagrangian (154) is symmetric under U(2) × U(1) c ; spinors transform in the fundamental representation of U(2), whereas U(1) c is associated with the chiral symmetry ψ → e iθγ5 ψ. The spectrum of the model involves a single SU(2) multiplet of massive fermions, with relativistic dispersion e(θ) = m cosh (πθ/2) and p(θ) = m sinh (πθ/2) (aside of the massless excitation charged under U(1) c that completely decouples). The amplitude of a two-fermion scattering is given by [96,164] Scattering of fermions carrying opposite spin is once again identical to the spin exchange of the SU(2) Heisenberg chain. Performing the particle-hole transformation on the 0th node assigned to physical (i.e. fermionic) density, i.e. Y 0 → 1/Y 0 , the resulting Dynkin TBA equations have the structure of the A 1 × A ∞ diagram, namely are of the form [164] log Y s = −δ s,0 β e + s log(1 Nodes with s ≥ 1 have been assigned to magnonic sstrings. Apart from an extra massive particle assigned to the initial node s = 0, the obtained equations are structurally strongly reminiscent of those of the isotropic Heisenberg spin-1/2 chain (see Fig. 9). A family of integrable Gross-Neveu models with SU(n) symmetry with n − 1 coupled copies of Eqs. (156) (A n−1 × A ∞ TBA incidence matrices) and Gross-Neveu models with O(2n) symmetry have been described in Ref. [150]. Numerical analysis. In order to extract the large-s asymptotic properties of thermodynamic state functions, we have numerically solved the TBA equations for the above SU(2)-invariant IQFTs with a common magnon structure. The deduced scaling properties match those of their spin-chain counterparts, as specified by Eqs. (125) and (126), indicating once again superdiffusive transport with dynamical exponent z = 3/2.

VII. CONCLUSION
We devoted this work to a systematic study of anomalous charge diffusion in equilibrium at finite temperature, found recently in certain quantum integrable systems with isotropic interactions. We aimed to settle, from a group-theoretical perspective, whether all integrable models exhibiting non-abelian continuous Lie group symmetries reveal the same type of transport anomaly. By establishing the link between the thermodynamic quasiparticle content and representations of the corresponding quantum group, we obtained universal algebraic dressing equations -a scaling analysis of these dressing equations, within the framework of generalized hydrodynamics, led to our conclusion that anomalous charge transport with a superdiffusive dynamical exponent z = 3/2 is generic to all integrable lattice ferromagnets described by Hamiltonians invariant under global non-abelian continuous symmetry, irrespective of the symmetry group or unitary irreducible representations associated with (local) physical degrees of freedom. Each such subclass indeed constitutes an infinitely large family of commuting integrable Hamiltonians, all of which exhibit z = 3/2 superdiffusion. We have argued that any other anomalous exponent is incompatible with the computed quasiparticle structure. In addition, the same type anomalous charge transport persists even in Lorentz invariant integrable quantum field theories with group-valued Noether currents. Due to its remarkable level of robustness, we dubbed this phenomenon as superuniversal.
Let us briefly restate the key steps that lead to this general conclusion. Recall, first, that anomalous transport occurs only when equilibrium density matrices possess full G-invariance: if one considers generic polarized states (realized by applying chemical potentials to the Cartan charges) the anomalous conductivity gets regularized and one recovers ballistic transport with subleading diffusive corrections, as one generically expects within generalized hydrodynamics. If one starts slightly away from the Gsymmetric state and approaches it, one finds that the Drude weight for ballistic transport vanishes, the diffusion constant diverges, and the quasiparticles that dominate the magnetic susceptibility and transport are macroscopically large coherent bound states of magnons, which we have called giant quasiparticles. Non-giant quasiparticles are effectively depolarized and do not contribute either to transport or to susceptibility. Transport of the Noether charges is in effect described by a dense gas of interacting giant quasiparticles. Each of these giant quasiparticles moves with a characteristic velocity inverse to its width: we established this result by explicit analysis of the GHD equations, and interpreted it as indicating that these giant quasiparticles are solitons made up of nonlinearly interacting, quadratically dispersing Goldstone modes. When these modes are present at finite density, they dress each other's charge in a nonperturbative way, and this nontrivial dressing leads to the fractional dynamical exponent z = 3/2.
To further elucidate the nature of the giant quasiparticles, we carefully examined the structure of semiclassical eigenstates from the viewpoint of an effective low-energy theory with respect to a (continuously-degenerate) ferromagnetic vacuum. At the level of an integrable classical field theory, giant quasiparticles can be identified with soliton modes, representing the nonlinear counterparts of quadratically dispersing (type-II) Goldstone bosons resolving linear fluctuations above the ferromagnetic vacuum. The number of internal (polarization) degrees of freedom of classical soliton modes is not equal to number of distinct flavour of the elementary excitations in the quasiparticle content of an integrable quantum chain as one would naïvely expect, but indeed exceeds the rank of the group G. This owes to formation of emergent classical multiflavored degrees of freedom called stacks produced by gluing together magnons of different flavors (which dissolve upon introducing quantum corrections). Stacks that contain momentum-carrying magnons should be regarded as independent physical excitations and we outlined how they can be naturally arranged on vertices of the Hasse diagrams. The total number of physical stacks is found to be in perfect agreement with the prediction of the non-relativistic variant of the Goldstone theorem. As a byproduct of our work, we provided a full classification of ferromagnets with the global symmetry of a simple Lie algebra and fundamental onsite degrees of freedom.
The anomalous nature of charge transport can be explained, at the most formal level, by invoking fusion identities amongst characters of Yangian symmetry algebras, which at thermodynamic quasiparticle spectra translates to a particular rigid algebraic structure. We thus wish to stress that the list of models discussed in this study is unlikely to be exhaustive and other classes of integrable models not included in our study are likewise expected to host a superdiffusive charge transport with characteristic exponent z = 3/2: we anticipate that the phenomenon extends beyond Yangians associated with simple Lie algebras to other classes of models based on rational solutions to the Yang-Baxter equations which includes (but is necessarily limited to) the Fermi-Hubbard model [36] and to other supersymmetric su(n|m) (or orthosymplectic osp(n|2m)) integrable spin chains that realize Yangian symmetries associated with Lie superalgebras (see e.g. [99,165,166]) whose bosonic Noether charges are known to display divergent diffusion constants [44], or integrable models based on infinite-dimensional affine symmetry algebras [167] (such as the Izergin-Korepin model [168]).
To benchmark our predictions, we have performed numerical simulations on a handful of representative instances of integrable quantum chain with local ferromagnetic exchange interaction, confining ourselves to classical simple Lie algebras and only to local degrees of freedom thattransform in the defining representation. We confirmed the predicted anomalous algebraic dynamical exponent z = 3/2 with great numerical precision. On the flip side, we did not manage to reliably discern the anticipated KPZ scaling profiles. We note that spin chains with symmetry of exceptional Lie algebras have been deliberately excluded from the present numerical study largely for technical reasons: on the one hand, algebraic construction of integrable quantum chains based on the 'quantized symmetry' of an exceptional Lie algebra is quite laborious (see e.g. [169]) and thus we find it better suited for a separate technical consideration. Another practical obstacle is that dimensions of the defining irreducible representations are (possibly with exception of (7) and (14) of G 2 ) conceivably too large to be efficiently simulated with current DMRG techniques.
Finally, a central piece of the full puzzle is still missing: despite integrability, a first-principles proof of the KPZ scaling profiles of dynamical structure factors remains elusive. At this time, there only exists a phenomenological picture which, under certain plausible assumptions, justifies the emergence of a noisy Burgers equation for the basic SU(2) case [43]. One obvious drawback of such an approach is that it cannot yield quantitative predictions such as the value of the KPZ coupling. Understanding quantitatively, from ab-initio principles, how KPZ uni-versality emerges (and not only z = 3/2) remains a major challenge for future works, even in the case of the SU(2) symmetric Heisenberg XXX spin chain. Another open problem is to appropriately extend the proposal of [42] to systems invariant under symmetries of higher-rank in order to extract the KPZ nonlinearity constant, which likely requires one to derive a semiclassical scaling limit for the NBA equations by fully accounting for presence of stack condensates.
for α ∈ ∆. Weights associated with the adjoint representation of g are identified with roots, that is E α are eigenvectors of the adjoint action ad H , [H i , E α ] = α(H i ) E α with α i ≡ α(H i ) denoting the ith component of a positive root vector α ∈ ∆ + . Roots are thus linear functional of t, and since every root is a linear combination of the simple ones α i , the latter provide the basis of t * dually paired to t.
Classification. Geometric structure of the root lattice is encoded in the r-dimensional Cartan matrix K, defined as with t a ≡ 2/(α a , α a ) ∈ {1, 2, 3}. The matrix elements of K are always integers. For simply-laced g we have t a = 1, whereas for non-simply laced algebra t a = 1 (t a = t ∈ {2, 3}) corresponds to long (short) simple roots. The Cartan matrices can be visualized by Dynkin diagrams (as depicted in e.g. Fig. 3): we draw an node (circle) for each diagonal elements K aa = 2, a = 1, . . . r; each two nodes are connected by 0 ≤ max(|K ab |, |K ba |) ≤ 3 lines (i.e. not connected when K ab = K ba = 0); when |K ab | > |K ba | we put an arrow pointing towards node b.
The corresponding highest-weight vector is annihilated by all the raising (Chevalley) generators E α for α ∈ ∆ + . Acting with the lowering generators E −α on a state with weight λ, its weight gets reduced by α.

Appendix B: Dressed charge fluctuations
Consider a general equilibrium macrostate characterized by quasiparticle densities ρ A (θ). Dressed charge susceptibilities within the Cartan sector can be linked to density fluctuations via Variations of q (i)dr A with respect to quasiparticle densities ρ A (θ) reads where K dr AB (θ) are the dressed scattering differential phases. Using that [174] δρ A (θ)δρ A (θ ) = 1 C AA (θ, θ ), and employing for convenience the following compact matrix notations where Ω[n] ≡ (1 + K n) −1 , we arrive at where q dr stands for the diagonal matrix of dressed charges q dr A . For the Cartan charges Q (i) considered in this paper, in an unpolarized background (i.e. h i = 0 for all i = 1, 2, . . . , r), we have q (i)dr a,s = 0, with exception in the limit s → ∞ where the above expression simplifies to give Eq. (25).

Appendix C: Kirillov-Reshetikhin modules
We consider an integrable quantum chain of length L with a nested spectrum. In order to find the number of distinct Bethe eigenstates in the finite-volume spectrum we consider an L-fold product representation of W-modules W Λp and its decomposition into g-modules V Λ (including non-trivial multiplicities) In the decomposition one finds multiplets with the highest weights which can occur with higher multiplicities. Here integers {N i } r i=1 denote the total number of Bethe roots of each flavor i ∈ I r , providing a bijective correspondence between the highest weight Λ and {N i } which determine a sector with fixed values of U(1) charges. Creating a quasi-particle of type a therefore reduces the maximal weight L Λ p of the reference ferromagnetic vacuum by the corresponding simple root α i .
Irreducible representations V Λ which appear in decomposition (C1) occur with multiplicities Z Λp ({N i }). These can be computed with aid of combinatorial formulae, first conjectured by Kirillov and Reshetikhin [80,82] (afterwards extensively investigated in a series of papers [85,175,176]), reading Here the sum runs over all integer sets (partitions) {n a,s } such that N a = s≥1 s n a,s , n a,s ≥ 0, 1 ≤ a ≤ r, and P a,s ≥ 0 for all a and s. Remark. While the described partitioning scheme assigns a unique pattern of quantum numbers to each eigenstate of the highest-weight type, it is remarkable that the is no bijective mapping with bound-state excitations (Bethe string). The correspondence is spoiled by exceptional solutions which involve so-called collapsed strings. The precise correspondence can nonetheless be established within the Wronskian formulation of Bethe Ansatz equations, as recently formulated in [111].

Multiplicities
The aim of this section is to illustrate the 'multiplicity formula' on two simple examples. We begin with A 1 ≡ su(2) of rank r = 1, with a single type of magnons in the spectrum. Within each fixed N -particle sector, each highest-weight eigenstate is characterized by a pattern of s-strings {n s }, such that N = s≥1 s n s . The multiplicity factor is where P s = L − 2 s ≥1 min(s, s )n s . For instance, in the fundamental SU(2) chain (the Heisenberg XXX model), Z ω1 (N ) for 0 ≤ N ≤ L/2 gives the numbers of distinct highest-weight Bethe eigenstates in the charge sector with N magnons that carries weight Λ = (L − 2N )ω 1 . The latter coincides with the branching coefficients in the decomposition of V ⊗L ω1 , that is b N = L N − L N −1 . As a concrete example we consider for the homogeneous spin-1/2 chain with L = 5 sites, whose Hilbert space decomposes into the following multiplets H ∼ = V ⊗5 ω1 = 5(2) + 4(4) + (6).

Classical characters
For the purpose of the next section, we collect the Weyl determinant formulae for characters of the classical simple Lie algebras g which can be found in standard text such as [172] or [173]. Knowing closed-form expressions for the classical g-characters proves helpful in finding explicit solutions to the constant T -system relations (are provided blow in section C 3) using the 'children expansion' to decompose the Kirillov-Resthetikhin modules over g-modules.
Constant T -functions encode grand-canonical Gibbs equilibrium state in the β → 0 limit. They thus depend solely on the U(1) chemical potentials x i associated with r Cartan charges of g. For simply-laced algebras, that is A r and D r and E r=6∼8 , the constant T -system assumes the form T 2 a,s = T a,s−1 T a,s+1 + b∼a T b,s , where ∼ stands for the neighboring nodes in the Dynkin diagram, that is the sum run over all b for which K ab = −1. The A-series is special in that W-modules are also irreducible as g-modules and therefore the T -functions are given by classical g-characters χ a,s themselves.
There are analogous formuale for the exceptional algebras which can be found in e.g. [83].

Character expansions
The T -functions are completely determined by the structure of the underlying classical Lie algebra g. Formally, the thermodynamic T -functions are identified as 'quantum characters' of the Kirillov-Resthetikhin modules W a,s associated to rectangular type representations of Yangians Y (g).
When viewed as g-modules, W a,s are generically not irreducible. The case of unitary algebras g = A n are exceptional since the W-modules W a,s are also irreducible as g-modules, implying that T -functions are precisely the su(2) characters χ a,s of rectangular irreducible representations V s ωa ≡ V a,s , cf. Eq. (106). For general g such a correspondence no longer holds as W a,s decompose nontrivially in terms of g-modules V a,s . It is known that such W-modules admit an expansion over g-characters χ(V λ ) ≡ χ λ of the form for some branching coefficients b λ . Decomposition (C28) is sometimes suggestively referred to as the 'children ex-pansion'. The general structure of such decomposition has been first conjectured in Refs. [81,82] and subsequently proved in a number of follow-up papers such as e.g. [176][177][178].
Using that functions T a,s satisfy the Hirota functional relation (110), the higher T -functions can be recursively constructed given the 'initial' ones T a,1 associated with characters of ath fundamental W-module W a,1 . For A r and C r≥2 , we have simply T a,1 = χ(ω a ). (C29) For B r≥2 and D r≥3 on the other have, we have T a≤r−1,1 = χ(ω a ) + χ(ω a−2 ) + . . . , (C30) T r,1 = χ(ω r ), and T a≤r−2,1 = χ(ω a ) + χ(ω a−2 ) + . . . , T r−1,1 = χ(ω r−1 ), T r,1 = χ(ω r ), respectively. There is a nifty graphical representation behind expansions of the form (C28) which can be carried out at the level of Young diagrams: for the C r ≡ sp(2r) case, the children diagrams are obtained from the parent rectangular diagram a × s by successively removing 1 × 2 'dominos', whereas for B r ≡ so(2r + 1) one likewise removes horizontal 2 × 1 dominos. One has to pay special attention to the presence of 'half-partitions' in the case of B-series, as we clarify below on an explicit example. Explicit combinatorial expressions for T a,s in terms of appropriately restricted sum can be found in Ref. [83]. In the unitary case Ar, the W-modules are already indecomposable. In the case of {Br and Dr}, they simply decompose in terms of g-modules with summands obtained by removal of 2 × 1 dominos (provided a < r in the case of Br, cf. remarks in appendix C 4). Similarly, for the C-series, one successively removes 1 × 2 dominos.

Examples
We conclude by providing a few explicit examples of the 'children expansion' of Yangian characters T a,s in terms of g-characters χ a,s .
As our first example we consider algebra B 2 = so(5) of rank 2. The initial (i.e. fundamental) T -functions read whereas the higher T -functions are given by the following restricted sums The first sum indeed involves a single term, yielding T 1,s = χ(s ω 1 ). The situation with the 2 × s quantum characters T 2,s is more subtle since the corresponding highest weights Λ = [0, s] are associated with so-called half-partitions Λ = ( s 2 , s 2 ) made of basic 1 × 1 2 rectangles. The children characters, which are obtained by successive removal of 2 × 1 dominos, therefore involve four basic rectangles. The consequence of this is a staggered odd-even structure.