Spin Drude weight for the integrable XXZ chain with arbitrary spin

Using generalized hydrodynamics (GHD), we exactly evaluate the finite-temperature spin Drude weight at zero magnetic field for the integrable XXZ chain with arbitrary spin and easy-plane anisotropy. First, we construct the fusion hierarchy of the quantum transfer matrices (T-functions) and derive functional relations (T- and Y-systems) satisfied by the T-functions and certain combinations of them (Y-functions). Through analytical arguments, the Y-system is reduced to a set of non-linear integral equations, equivalent to the thermodynamic Bethe ansatz (TBA) equations. Then, employing GHD, we calculate the spin Drude weight at arbitrary finite temperatures. As a result, a characteristic fractal-like structure of the Drude weight is observed at arbitrary spin, similar to the spin-1/2 case. In our approach, the solutions to the TBA equations (i.e. the Y-functions) can be explicitly written in terms of the T-functions, thus allowing for a systematic calculation of the high-temperature limit of the Drude weight.


Introduction
Quantum integrable systems are characterized by an infinite number of local or quasi-local conserved charges, which have led to extensive investigations into their unique physical properties.While conventional systems are understood through Gibbs ensembles, integrable systems, due to their multitude of extensive conserved charges, are described by Generalized Gibbs Ensembles (GGE) [1,2] (refer to [3] for a review).This distinction becomes particularly evident in thermalization processes, where integrable systems often deviate from the predictions of frameworks like the eigenstate thermalization hypothesis [4].
On the other hand, our understanding of non-equilibrium properties in integrable systems has seen rapid advancements, largely attributed to the theoretical framework of Generalized Hydrodynamics (GHD) [5,6] (see also [7] for lecture notes).GHD effectively merges the concepts of fluid mechanics and GGE, providing a novel perspective in studying the nonequilibrium properties of quantum integrable systems (recent reviews are available in [8][9][10]).Moreover, the GHD framework has been extended to weakly broken integrable systems (refer to [11,12] for a review).
Among these advancements, the finite-temperature spin transport properties of the spin-1/2 XXZ chain, which are deeply relevant to this work, have been extensively studied for a considerable period [13].The spin-1/2 XXZ chain is given by where σ x j , σ y j , and σ z j are the spin-1/2 spin operators (Pauli matrices), and ∆ stands for the anisotropy parameter.In this model, the spin current is not a conserved quantity except for ∆ = 0, which corresponds to a free fermion system.For |∆| < 1, the system exhibits ballistic spin transport, characterized by a finite spin Drude weight [14][15][16][17] and also by a dynamical exponent z = 1 that signifies a scaling relation x ∼ t1/z = t between the magnetization displacement x and time t.In contrast, when ∆ > 1, the spin transport is believed to be diffusive with a dynamical exponent z = 2 [18,19].Interestingly, at the isotropic point ∆ = 1, i.e., at the boundary between the diffusive and ballistic regimes, the spin transport at infinite temperature is considered to exhibit superdiffusive behavior with a dynamical exponent z = 3/2 [18,20,21].Further investigations [22] suggest that, particularly at this isotropic point, the transport properties are likely to belong to the Kardar-Parisi-Zhang (KPZ) universality class [23].On the other hand, a recent investigation [24] indicates that the transport at ∆ = 1 cannot be explained by the KPZ type, which has also been supported by recent quantum simulations [25].From this perspective, the topic continues to be a subject of ongoing debate.More rigorous and quantitative analysis of spin transport at the isotropic point is desired, especially a more accurate evaluation of the exponent.
For |∆ = cos(π/p 0 )| < 1, the spin Drude weight at finite temperatures had long been a topic of debate since Zotos derived it in [14] by combining Kohn's formula [26] with the Thermodynamic Bethe Ansatz (TBA) [27].Taking a notably different approach, Prosen and Ilievski evaluated the Drude weight at infinite temperature [15].Specifically, they constructed quasi-local conserved charges and optimized the Mazur bound of the Drude weight using them.Remarkably, this optimized lower bound, referred to as the Prosen-Ilievski bound, exhibits a fractal-like structure, more aptly described as popcorn-like [28] 1 .Moreover, with the emergence of GHD, the spin Drude weight has been formulated in a universal manner in [29,30].The result from GHD reproduces Zotos's formula.Intriguingly, the high-temperature limit of the Drude weight exactly coincides with the Prosen-Ilievski bound [16].This popcorn structure is intrinsically linked to the structures of strings [27] that define the thermodynamics of the XXZ model.In fact, recent studies suggest that similar popcorn structures also appear in higher-spin integrable XXZ models and the sine-Gordon model [28,31,32], where these thermodynamics are governed by the same strings as in the XXZ model.
This paper focuses on the finite-temperature spin Drude weight at zero magnetic field for an integrable XXZ chain with arbitrary spin S. We investigate the model in the critical regime, characterized by the anisotropy parameter ∆ = cos(π/p 0 ) where p 0 ∈ Q ≥2 .This model was introduced by Kirillov and Reshetikhin [33] through a "fusion" of the XXZ model.For S = 1, the model corresponds to the Zamolodchikov-Fateev model [34].Regarding the spin Drude weight, a suboptimal Mazur bound was derived for the high-temperature limit in the S = 1 case [35].More recently, the exact high-temperature behavior of this quantity has been further clarified for the specific case where S = 1 and p 0 ∈ N ≥2 [28].
As reported in [33,36], for S ≥ 1, the anisotropy ∆ must be constrained by the value of S to ensure that the Hamiltonian derived from the fusion procedure possesses real eigenvalues.Specifically, for an integer spin S greater than 2, the permissible region of ∆ is composed of several separate open intervals.This also holds for half odd integer spin S greater than 7/2.(See Table 1 in the main text.)Importantly, the low-energy properties in each region, for a given S, are governed by a distinct conformal field theory [36,37].From this perspective, it is crucial to delve into how this rich phase structure affects the transport properties.
We investigate the spin Drude weight of this model using GHD across a broad range of temperatures and anisotropies.To systematically derive the exact expressions at infinite temperature, we derive the TBA equations using the quantum transfer matrices (QTMs) and their functional relations referred to as the T and Y -systems [38].Notably, the resulting Drude weight, at both finite and infinite temperatures, prominently displays a popcorn structure that depends on ∆ and exhibits discontinuities throughout.As expected from the different excitation properties in each region explained earlier, the behavior of the Drude weight also varies distinctly across regions.
The paper is organized as follows.In the subsequent section, we provide a comprehensive overview of the integrable XXZ chain with arbitrary spin, including its derivation using the fusion procedure.In section 3, we construct the QTMs and derive their functional relations, called the T and Y -systems.By analyzing their analytic structure, we derive the TBA equations.In section 4, using the solutions to the TBA equations, we calculate the finitetemperature Drude weight at zero magnetic field over a wide range of anisotropies and spins.Section 5 presents the exact expression for the Drude weight in the high-temperature limit.Section 6 is devoted to a summary and discussion.Technical details and supplementary information necessary for the main text are provided in the Appendices.

Integrable XXZ chain with arbitrary spin
In this section, we construct the integrable XXZ chain with arbitrary spin S = σ/2 (σ ∈ N + ), applying the so-called fusion procedure to the six-vertex model.The energy spectrum of the model is derived by means of the Bethe ansatz as in [33,36].

Higher spin generalization of the six-vertex model
The six-vertex model is a classical statistical model defined on a two-dimensional square lattice, where the configuration variables (classical spins) {1, 2} are assigned on the edges between each vertex of the lattice.Admissible configurations of the model are described by the following local states with the Boltzmann weights where (2. 2) The weights for all other configurations are set to zero.Let V be a two-dimensional vector space spanned by an orthonormal basis {|1⟩, |2⟩}.
satisfies the Yang-Baxter equation (YBE): Here R jk (v) stands for the operator R(v) acting non-trivially on where V j is an indexed copy of V .Graphically, the YBE (2.4) is depicted as This model is related to the spin-1/2 XXZ chain as seen later ((2.22) and (2.28)).
Starting from the six-vertex model, we can construct a higher spin version for the sixvertex model by the fusion procedure [39] (see also [40] and [41]).Let where σ − j is the lowering operator acting on the jth space V j as σ (2.9) One finds that In fact, for the case that |w⟩ ⊗ |x⟩ ∈ V (m) ⊗ V (n) and R jj+1 acts on |x⟩, we find (2.10) for all 1 ≤ j ≤ n−1, where P jj+1 is the permutation operator on V j ⊗V j+1 .Note that the first equality in the above follows from the YBE (2.6), and the second equality is the consequence of |x⟩ ∈ V (n) .Eq. (2.10) Thus, we define the fusion R-matrix as R (m,n) (v) restricted to the subspace V (m) ⊗ V (n) , and denote it by The elements of R where α, γ ∈ {1, 2, . . ., m + 1} and β, δ ∈ {1, 2, . . ., n + 1}.Explicitly, they read (2.13) In the above, the position of the δ − 1 (resp.γ − 1) 2's on uppermost (resp.rightmost) edges can be taken arbitrarily, since V (m) and V (n) are preserved under the action of R (m,n) (v) as explained before.By repeatedly applying the YBE (2.6), one can show that the fusion R-matrix also satisfies the YBE: Here, let us provide some remarks about the fusion R-matrix R (m,n) (v) constructed by (2.13).The R (m,n) (v) possesses the following symmetries (C and P invariances): ( The charge conservation also holds: In general, our construction does not guarantee T invariance.Namely, for certain combinations of α, β, γ, and δ.This is in contrast to the R-matrix used in [33,36] (see also [42,43]), where the fusion procedure is designed to ensure the T invariance of the R-matrix.However, a similarity transformation renders a T -invariant matrix [33,44,45], which coincides with the R-matrix in [33,36].Here, is not necessary in the following arguments, we list them for a few small values of m: Since the similarity transformation does not change the eigenvalues of the transfer matrices, in this paper, we will develop the argument based on the R-matrix (2.13) instead of (2.18), whose elements can be intuitively evaluated graphically.

Integrable XXZ chain with arbitrary spin
The integrable XXZ model with arbitrary spin S = σ/2 (σ ∈ N + ) can be derived by first taking the logarithmic derivative of the row transfer matrix of the fusion six-vertex model with respect to the spectral parameter v ∈ C and then applying a suitable transformation to the resulting expression to make it a Hermitian operator.Let T R n (v) ∈ End V (σ) ⊗L be the row transfer matrix: As a consequence of the YBE (2.14), the transfer matrices commute for different values of the spectral parameter and different fusion levels Let us define the Hamiltonian H (σ) ∈ End V (σ) ⊗L for the spin-σ/2 integrable XXZ model, defined on a 1D periodic lattice with L sites (i.e., V where jj+1 is the permutation operator on j+1 , σ z j is the z-component of the spinσ/2 operator acting on the state |σ; J ≶ 0 is a coupling constant, and h is the external magnetic field.The second equality in (2.22) comes from which is a consequence of R jk (0 due to the charge conservation (2.16).It is crucial to point out that the operator ‹ H (σ) is, in general, non-Hermitian (more precisely, a real asymmetric matrix) except for σ = 1 due to the breaking of the T invariance for the R-matrix (see (2.17)).However, a similarity transformation corresponding to (2.18) reproduces the Hamiltonian H (σ) given by [33,36]: where is a diagonal matrix.In particular, for S = 1/2 (σ = 1), ‹ H (1) = H (1) , which is Hermitian by definition of the six-vertex model, is the well-known spin-1/2 XXZ model: (2.28) The Hamiltonian for the higher-spin integrable XXZ chain is not simply obtained by replacing the spin operators σ x,y,z j in (2.28) with those of a higher spin representation.(Note that [σ a j , σ b j ] = 2iϵ abc σ c j , where a, b, c ∈ {x, y, z}, and ϵ xyz = 1 is a completely antisymmetric tensor.)For instance, for S = 1 (σ = 2), we can derive ‹ H (2) where σ x,y,z j are the spin-1 operators, σ j := (σ x j , σ y j , σ z j ) and (2.30) The operator ‹ H (2) 0 defined above is obviously non-Hermitian, however it can be transformed to the Hermitian operator H (2) 0 by the similarity transformation (2.26), where F (2) j is explicitly given by (2.19).The resultant H (2) 0 is obtained by replacing in (2.29), which exactly coincides with the Zamolodchikov-Fateev model [34].Up to authors' knowledge, the explicit form of the Hamiltonian has been found only for the cases of σ = 1, σ = 2 [34] and σ = 3 [46].The eigenvalues of T R n (v) (2.20) and ‹ H (σ) are preserved under the similarity transformation (2.18) and (2.26).Also, the expectation values and the matrix elements of any operator o with respect to the eigenstates of H (σ) coincide with those for the corresponding operator with respect to the eigenstates of ‹ H (σ) .For instance, the equilibrium thermal averages Tr(e −β ‹ H (σ) ) , (2.33) are exactly the same, i.e., ⟨o⟩ = ⟨ o⟩. (2.34) Physically, the spin Drude weight discussed in this paper is the higher-spin integrable XXZ chain (2.26), which can be evaluated by using the model (2.22) through the correspondence as in (2.34).
The eigenvalues of the row transfer matrix T R n (v) can be obtained by the analytic Bethe ansatz [47], which is a shortcut procedure that allows one to find the eigenvalues of a transfer matrix without knowing the details of the eigenstates.The resultant eigenvalues (denoted by the same notation T R n (v)) is written in the dressed vacuum form (DVF) [33]: where Note that M = 0 (i.e., q(v) = 1) corresponds to the vacuum state which is an eigenstate of T R n (v).The vacuum eigenvalue is actually calculated by the substitution of which easily follows from (2.13), into q ⟨0|T R n (v)|0⟩ q .The M -particle states are spanned by the basis set The unknown {λ j } is determined by the Bethe ansatz equation (BAE): Up to the overall factor and the constant term (the second term), the energy spectra agree with those obtained in [33,36].
In this paper, we consider the case where the model is in the critical regime, i.e., In this case, it should be noted that for a given σ, except for σ = 1, the anisotropy parameter ∆ cannot necessarily take any value in the range 0 ≤ ∆ < 1. Namely, for a given σ, only those values of the parameter θ that satisfy the condition v(σ + 1) sin(θj) sin (θ(σ + 1 − j)) > 0 for all j ∈ {1, 2, . . ., σ}. (2.43) Here v(σ + 1) is the spin parity, where v(n) is defined as and ⌊x⌋ denotes the largest integer less than or equal to x. See Table 1 for the admissible region of the parameter p 0 for several given σ.
Originally, the relationship (2.43) between the spin σ/2 and the anisotropy parameter ∆ = cos θ was found by Kirillov and Reshetikhin in [48,49].They discovered this by considering the conditions for the existence of string-type solutions to the BAE (2.40), which are given by the form Here, δ > 0, n j ∈ N + is the admissible string length known as the Takahashi-Suzuki (TS) number defined by (A.14), and is the center of the string.The integers v j ∈ {1, −1} and v(σ + 1) ∈ {1, −1} are, respectively, the string parity defined as (A.17) and the spin parity (2.44).These can also be rewritten using the sequences w j and w j defined as (A.15) as v j = (−1) w j and v(σ + 1) = (−1) ‹ w jσ +1 , where j σ is later defined as in (3.23).The condition (2.43) is equivalent to the requirement that the admissible string lengths {n j } must contain σ + 1: (To be more precise, this condition is valid when p 0 is an irrational number.On the other hand, when p 0 is a rational number, the condition becomes (3.22).See also Appendix B for details.) Subsequently, Frahm, Yu, and Fowler found that the conditions in (2.43) are simply those necessary for the Hamiltonian to possess real eigenvalues [36].Clearly, for σ = 1, the arbitrary θ ∈ (0, π/2] satisfies (2.43).On the other hand, for a given σ ≥ 2, (2.43) is fulfilled by some open intervals for θ (see Table 1 for an example and [36] for details).
The ground state, low-lying excitations, and finite-temperature properties of the model for arbitrary (σ, θ) satisfying (2.43) have been investigated in [33,50].This was done by analyzing the thermodynamic Bethe ansatz (TBA) equations derived in [33].Specifically, f , representing the free energy per site, is given by where p j (v) is the bare momentum, explicitly defined in [33,36].Additionally, is a sign factor, with q j defined in (A.18).The function η j (v) is determined by the TBA equations log where f * g(v) represents the convolution defined by dx, and ε j (v) is the bare energy defined as For the explicit form of the integration kernel K jk (v) in (2.49), see [33,36].The conformal properties have also been studied in [36,37,51,52].Characteristically, for some intervals (σ, θ) satisfying (2.43), the ground state is described by multiple types of string solutions.Thus, there are multiple branches of low-lying excitations with different Fermi velocities within these intervals, and as a result, the free energy per site f at low temperatures T behaves like where e 0 is the ground state energy per site, v f j is the Fermi velocity for the n j -string describing the ground state, and c j is the central charge, which is, in general, given by the form of c j = 3k j /(k j + 2) (k j ∈ N + ).The integer k j is determined through a rather complicated procedure in [50].For example, for p 0 in one of the intervals k < p , the ground state in this interval is described by a single type of string, and then the central charge is given by c = 3k/(k + 2) [36].

QTM and TBA
Let us formulate the model at finite temperatures.The bulk thermodynamic quantities are described by the TBA equations, a set of coupled non-linear integral equations.Concerning transport properties, the Drude weights are universally given by using the solutions to the TBA equations [29,30].The TBA equations for the present model have already been established in [33], whose explicit form is given by (2.49).In this section, however, we revisit the equations using the quantum transfer matrix (QTM) approach [53][54][55][56][57][58][59][60][61] (see also books [62][63][64]), enabling us to develop a systematic analytical formulation for the Drude weight at high temperatures.

QTM
To consider the thermodynamics of the original quantum system, let us go back to the twodimensional classical system, i.e., the higher-spin generalization of the six-vertex model whose local Boltzmann weights are defined as (2.13).On the N × L square lattice, we define the following double-row transfer matrix where the matrix depicted by the dashed line (we denote it by R Using (2.22), (2.24), (3.2) and (3.3), T DR σ (u, 0) can be expanded with respect to u as which gives the partition function Z.That is, by setting the free energy per site f can be formally expressed as Here N is the Trotter number N ∈ 2N + , and β is the reciprocal temperature: β := 1/T .To actually evaluate (3.6), we now introduce the QTM T n (u, v) ∈ End V (σ) ⊗N : See also Fig. 1.Note that the operator σ z 1 in (3.7) is the z-component of the spin-n/2 operator acting on (see (2.23) for the definition of σ z j ).The free energy (3.6) is then given by the largest eigenvalue of T σ (u N , 0), denoted by Λ(u N , 0), as Here we have used the exchangeability of the two limits in the first equality as proved in [54,55] and the fact that there exists a finite gap between the largest eigenvalue Λ(u N , 0) and the subleading eigenvalues.Applying the YBE (2.14), we can show that the QTMs commute with each other for different values of v and different levels of n, as long as the parameter u N is the same: (3.9) The eigenvalues of T n (u N , v), also denoted by the same symbol T n (u N , v), are written in the DVF by the analytic Bethe ansatz.From (3.9) and the fact that the algebraic structure of ) is identical to that of the fusion QTM for σ = 1 defined in [38], the DVF is the same as for σ = 1 except for the vacuum part.The vacuum state |0⟩ a of T n (u N , v) is given by (2.37) for the row transfer matrix), where the state |m; n⟩ is defined by (2.7).On the other hand, the M -particle eigenstates of the QTM are spanned by the basis set Figure 1: The QTM T σ (u, v) and the double-row transfer matrix T DR σ (u, v) defined on the N × L square lattice, where N ∈ 2N + is the Trotter number.Noticing (2.23), and substituting the following elements of fusion R-matrix which can easily be derived from (2.13), one evaluates the vacuum eigenvalue a ⟨0|T n (u, v)|0⟩ a .The resulting DVF (henceforth referred to as the T -function) is given by Note that the dress part in the above formula is also equivalent to that for the row transfer matrix (2.35) after changing the variables as v → −v, and v j → −ω j .The M unknown numbers ω j (1 ≤ j ≤ M ) are determined by the BAE: In the critical regime (2.42) restricted by the condition (2.43), the largest eigenvalue Λ(u N , 0) that gives the free energy of the model (see (3.8)) lies in the sector M = N σ/2.More explicitly, the BAE roots {ω j } M j=1 associated with the largest eigenvalue Λ(u N , 0) of the QTM T σ (u N , 0) are given in the following form: Here, ϵ k l ∈ R (1 ≤ l ≤ N/2) are symmetrically distributed about the imaginary axis for a given k, and δ k l ∈ R. In particular, in the high-temperature limit β → 0 (i.e., u N → 0), we observe that the BAE roots (3.15) are reduced to ω k l = i(σ + 1 − 2k).In fact, substituting them into (3.13)yields log Λ(u and hence Λ(u N = 0, 0) (3.17) Thus, from (3.8), − lim β→0 βf (i.e, the entropy per site at β = 0) reproduces the physically relevant result − lim β→0 βf = log(σ + 1).

TBA equation via T and Y -functions
Our first goal is to relate the solutions to the TBA equations in terms of the T -functions (3.13).As shown in (3.8), the thermodynamic quantities of the spin-σ/2 XXZ chain are described by the largest eigenvalue of the QTM T σ (u N , v), which is given by (3.13) via the solutions to the BAE (3.14).However, in general, the Trotter limit N → ∞ in (3.8) cannot be taken analytically by directly solving the BAE, except for the high-temperature limit T = ∞ (see (3.17)).To overcome this, we embed T σ (u N , v) into the functional relations (T and Y -systems) that are fulfilled by the T -functions and their suitable combinations called Y -functions.By analyzing analytic properties, the Y -systems are mapped to the form of nonlinear integral equations, in which the Trotter limit can be taken analytically.These nonlinear integral equations are the TBA equations.
In what follows, the common variable u in the T -functions is often omitted to simplify the notation.As previously described, the algebraic structure of V (n) of the QTM (3.7) is identical to that of the QTM for the spin-1/2 model.Consequently, the T -systems defined by (3.13) are the same as those obtained for the spin-1/2 model in [38]: which holds for any v ∈ C, for any ω j ∈ C (referred to as the off-shell roots, implying that there is no requirement for ω j to satisfy the BAE (3.14)), and for all integers n, y ∈ N + such that n ≥ y ≥ 1.Note here that and the T -functions have the periodicity: Now let us restrict ourselves to the case that p 0 defined in (2.42) (see also (2.43)) is a rational number and express it as the continued fraction form: where α ≥ 1, ν 2 , . . ., ν α−1 ∈ N + and ν α ∈ N ≥2 .We also set ν 1 ∈ N ≥2 which comes from the condition p 0 ≥ 2 (see (2.42)).Note that irrational p 0 corresponds to α = ∞.In Appendix A, we list the sequences required in this section: {m r } α r=1 , {p r } α r=1 , {y r } α r=−1 , {z r } α r=−1 , the TS-numbers {n j } j≥1 [27] defining the admissible lengths of the strings (2.45), and its slight modification { n j } j≥1 [38].We also introduce the sequence {w j } j≥1 and { w j } j≥1 in (A.15).related to the string parity (A.17).All of these sequences are uniquely determined for a given value of p 0 .In Appendix B, we prove that the restriction between the spin and the anisotropy, given by (2.43), is equivalent to the condition Hereafter, we consider the physical case where σ and p 0 satisfies (3.22).Correspondingly, let us define the index j σ ∈ {1, 2, . . ., m α − 1} and r σ ∈ {1, 2, . . ., α} such that Note that j σ=1 = 1 and r σ=1 = 1 for any p 0 satisfying (3.21).
With the definitions provided above, we introduce a functional relation involving T n (v), as described in [38]: where M denotes the number of BAE roots for (3.14).Note that this relation holds true for any off-shell roots, just as in equation (3.24).The proof can be straightforwardly completed using equation (A.8) and the following periodicities: It is important to note that the relation (3.24) is valid only when p 0 is an arbitrary rational number as given by (3.21).
The Y -functions {Y j (v)} mα j=1 are defined by the following combinations of the T -functions.For m r ≤ j < m r+1 (0 ≤ r ≤ α − 1) and j ̸ = m α − 1, while for j = m α − 1 and j = m α , We also set Y 0 (v) := 0 and Y −1 (v) = ∞.From the T -systems (3.18) and (3.24), we have for m r ≤ j < m r+1 (0 ≤ r ≤ α − 1) and j ̸ = m α − 1, and One finds that the above definitions of the Y -functions essentially the same with those introduced in [38] for the case where h = 0 and σ = 1.The sole distinction arises from a uniform shift in the argument that depends on σ, namely, i w jσ+1 p 0 , which reflects the spin parity for the strings (2.45).Consequently, the Y -functions satisfy the same Y -systems as given by [38]: the last two equalities in (3.30) reduce to Then, by newly introducing in (3.30), we reproduce the Y -system originally introduced in [38].In this paper, however, we use (3.30), because the physical structure of the quasi-particles, which is essential in the GHD formalism, is more clear.By substituting a set of the BAE roots (on-shell roots) that gives the largest eigenvalue Λ(u N , 0) of the QTM T σ (u N , 0), and examining the analytical properties of the Y -functions, the Y -systems are transformed into the following non-linear integral equations through the Fourier transform as explained in Appendix C.
In the non-linear integral equations (NLIEs), the Trotter limit N → ∞ can be taken analytically.The resultant NLIEs read where and with {q j } and { q j } defined by (A.18).We have numerically confirmed for wide ranges of σ and p 0 that the solutions η j (v) to the above TBA equations agree with those obtained by the string hypothesis [33] as given in (2.49).The free energy (3.8) can be expressed in terms of η j (v) as (D.12).See Appendix D for detailed derivation of the free energy.

Drude weight
Now we evaluate the finite-temperature spin Drude weight at zero magnetic field (i.e., h = 0).A general formula describing the Drude weight in terms of quantities determined by the generalized TBA equations (E.27) is explained in Appendix E, which is based on the arguments in [7,30].
In general, the spin Drude weight is defined by the Kubo formula using the dynamical correlation functions in the equilibrium state: where j s represents the spin current density, and ⟨•⟩ c denotes the thermal average of the connected correlation function: ⟨o 1 o 2 ⟩ c := ⟨o 1 o 2 ⟩ − ⟨o 1 ⟩⟨o 2 ⟩ with ⟨•⟩ defined as (2.33).Note that the quantity D s (β)/β is sometimes used as an alternative definition of the Drude weight.
In such cases, the factor β in front of the integral in (4.1) becomes unnecessary.However, in this paper, we adhere to the definition (4.1) as used in [14,16].For a generic expression of the Drude weight, see (E.1).The Drude weight D 00 (β) in Appendix E is denoted as D s (β) in this section.By projecting the currents onto the space spanned by the local or quasi-local conserved charges, the dynamical current-current correlation function, defining the Drude weight, can be expressed in terms of equal-time current-charge and charge-charge correlation functions [7,30] as given by (E.16).For the present system (2.26), the local spin current (j s ) kk+1 (1 ≤ k ≤ L) is determined by which follows from the discrete form of the continuity equation (cf.(E.2)): where ∂ t (q s ) k = i î H (σ) , (q s ) k ó .The corresponding quantities ( q s ) k and ( j s ) k (see (2.32)) are also given by where the second equality in the first equation is due to the fact that F (σ) and σ z k are diagonal matrices.
The averages of (4.2) over the GGE (see (E.6)) can be expressed in terms of a mode decomposition using the n j -string (1 ≤ j ≤ m α ) as shown in (E.17), (E.32), and (E.33).They reduce to where must hold, since ⟨•⟩ and ⟨ •⟩ defined as (2.33) denote the equilibrium thermal average (cf.(E.6)).In (4.5), the first equalities are derived from (2.34), ρ j (v) is the distribution function associated with the n j -string, and h s;j (v) (corresponding to h 0;j (v) in Appendix E) denotes the bare spin, which is the number of magnons composing the n j -string, i.e., h s;j (v) = h 0;j (v) = n j .Accordingly, by (2.49) and (E.27), we identify in the present case.The effective velocity, v eff j (v), is generally described by (E.31) via the generalized TBA equations (E.27).In the context of our present case, it is defined as where η j (v) is determined by the standard TBA equations (3.34), which are equivalent to (2.49).As v eff j (v) in (4.8) is an odd function while ρ j (v) and h s;j (v) are even functions, we can easily confirm (4.6).Then, using the procedure described in [7,30], which is also summarized in Appendix E, we can express the spin Drude weight as Here, h dr 0;j (v) represents the dressed spin, which is defined as as indicated in (4.7) and (E.29).Moreover, according to the TBA equations (3.34), for zero magnetic field (h = 0), the value of h dr 0;j (v) is , where ρ h j (v) is the hole density of the n j -string, into the fermionic distribution function ϑ j (v) (eq.(E.20)) with the density of states ρ tot j (v) := ρ j (v) + ρ h j (v), we obtain Here ξ j is a sign factor defined in (2.48).To derive the last equality, we have used (E.30) with ε dr j (v) = ∂ β 1 log η j (v) = ∂ β log η j (v) which follows from (E.29), (E.18) and (E.19).By substituting (4.8), (4.11) and (4.12) into (4.9), and taking into account along with the identity η mα (v) −1 = η mα−1 (v) which holds for h = 0, one sees that only the last two η-functions η α−1 and η α contribute to the Drude weight at h = 0. Explicitly it reads The formula presented above (4.14) is formally the same as that for the spin-1/2 model [14,16,17,65].The spin dependency is given through the η-functions (3.34).We solved the TBA equations numerically to compute the spin Drude weight (4.14).In Fig. 2, the temperature-dependence of the Drude weight (divided by σ 2 ) D s (β)/σ 2 is shown for the cases where the model with σ = 1, 2, 3, 4 (S = σ/2) in the antiferromagnetic regime J = 1.Although the dependencies for σ = 1 have already been known in [14,16], we included them for the sake of comparison.As mentioned earlier, the physically admissible region of ∆ is restricted as in (2.43) (refer to Table 1 for concrete examples).For a given ∆ and σ, the Drude weight varies smoothly with the temperature, and approaches zero at the isotropic limit p 0 → ∞ (∆ → 1).Notably, when the admissible region of ∆ consists of separate intervals (as with σ = 4 in the figure, where the admissible region of p 0 is separated into two intervals: p 0 ∈ (2, 3) and p 0 ∈ (4, ∞)), the behavior of the Drude weight differs ) and various anisotropies ∆ = cos(π/p 0 ).D s /σ 2 smoothly varies with respect to T and approaches zero at the isotropic limit p 0 → ∞ (∆ → 1).For σ = 4, the admissible region of p 0 is separated into two intervals: p 0 ∈ (2, 3) and p 0 ∈ (4, ∞) as in Table 1.The behavior of D s /σ 2 in each region distinctly differs, reflecting the fact that the structure of the energy spectrum characterized by the strings is essentially different in each interval.1 for details.
distinctly in each region.As pointed out in [33,36,37,50], the ground state and the lowlying excitations in each region are characterized by different types of strings.Accordingly, the low-energy properties in each region are described by a different conformal field theory [37].The different behavior of the Drude weight can be interpreted by these differences in the structure of the energy spectrum.In Fig. 3, we also depict the anisotropy dependence of the Drude weight for J = 1 divided by β = 1/T and σ 2 , i.e., D s /(β 2 σ 2 ).Again, for comparison, we include the known results [16] for σ = 1.For a given T and σ, the anisotropy dependence of the Drude weight exhibits the characteristic popcorn structure.The popcorn structure is suppressed with the decrease in temperature and will eventually exhibit a smooth curve depending on ∆ at the zero-temperature limit.

High-temperature limit
We aim to derive an analytical expression for the high-temperature behavior of the Drude weight at h = 0 (4.14), namely lim β→0 D s (β)/β.Given the complexity of the TBA equations, we take a different approach by adopting and extending the method developed in [16] for the spin-1/2 case.Namely, we evaluate this limit by directly examining the high-temperature behavior of the T and Y -functions, as defined in equations (3.13), (3.26) and (3.27).
In fact, according to the definition of free energy given by equations (3.6) through (3.4) and (3.5), considering only the case N = 2 would suffice for evaluating contributions of O(β).However, in this work, we extend the discussion to arbitrary even values of N to make the argument more generally applicable.The behavior of T n (u N , v) at O(β) can be evaluated by inserting the BAE roots, specifically those given by (3.15), up to terms of O(β).Since the real parts of the roots are symmetrically distributed with respect to the imaginary axis, we have where and z 0 l := −iu N and z σ+1 l Second, due to (5.1), we can rewrite , where the overall constant has been determined by the leading coefficient of z N/2 in (5.2).Then, by comparing the coefficients of z N/2−1 , we obtain Here we have used the equality N/2 l=1 ϵ k l = 0. Finally, imposing the condition δ 1 = −δ σ , which comes from the symmetrical argument, we arrive at where we introduce to simplify subsequent notations.Inserting equation (3.13), decomposing the summation appearing in T yα+y α−1 −1 (v ′ ) into two parts as , and using the periodicities (3.25), one obtains . (5.6) Due to the BAE (3.14) and periodicities given in (3.25), f (v) is an entire function and is bounded.Therefore, by Liouville's theorem, f (v) must be a constant which is calculated by taking the limit v → ∞.Using M j=1 ω j = 0, it explicitly reads (5.7) Taking into account (3.5) and (5.3), and expanding g(v) up to O(β), we obtain (5.8) Substituting (5.5) in the above and using where the function a j (v) is defined as sin(θq j ) ch(θv) + cos(θq j ) . (5.11) Finally, the insertion of (5.10) into (4.14)yields dv. (5.12) For S = 1/2 (σ = 1), (5.12) coincides with the result obtained in [16] up to an overall factor, which is due to the difference in the definition of the coupling constant J.In this case, the integral can be explicitly computed in [16], yielding: This result is consistent with the Prosen-Ilievski bound [15] derived via the construction of quasi-local conserved charges [66,67].On the other hand, for spins higher than 1/2, the integrand becomes more complicated.It appears difficult to calculate it explicitly when p 0 is a generic rational number.For S = 1 (σ = 2), the integral can be explicitly calculated for α = 1, i.e., p 0 = π/ν 1 .The result is which exactly coincides with the result recently derived in [28].In Fig. 4, we depict the high-temperature limit β → 0 of D s (β)/(σ 2 β).The result for σ = 1 is referred to as the Prosen-Ilievski bound [15].For σ ≥ 2, the spin Drude weights exhibit a prominent popcorn structure similar to the σ = 1 case.For p 0 close to, but not exactly on, the irrational value, which can be obtained with the desired precision by taking the limit α → ∞, the Drude weights are approximately represented by the points on the curve.Namely, the Drude weight behaves like a continuous function for p 0 close to the irrational value, but is discontinuous everywhere for rational p 0 .The Drude weight cannot be obtained in the present approach when p 0 is exactly an irrational number.

Summary and Discussion
We have examined the finite-temperature spin Drude weight for the integrable XXZ chain with arbitrary spin S at zero magnetic field.In the critical regime, the admissible string lengths are described by the TS-numbers, akin to the spin-1/2 case.In relation to this, the Drude weight exhibits a behavior similar to that observed in the spin-1/2 XXZ, displaying a pronounced popcorn structure that is discontinuous everywhere with respect to the anisotropy parameter ∆.For a given S, the admissible region of ∆, in general, comprises several separate intervals, which is a fundamental distinction from the spin-1/2 case.Importantly, the structure of the TBA equations in each region is distinct, resulting in the Drude for the anisotropy dependence.The graph displays popcorn structures, indicating everywhere discontinuity with respect to ∆.For the physically admissible region of ∆, refer to Eq. (2.43) and Table 1.
weight behaving drastically differently in each region.In this work, we have constructed the TBA equations, using the QTMs and their functional relations.As a significant advantage, this approach enables us to perform a systematic exact calculation of the spin Drude weight in the high-temperature limit.Let us briefly discuss the open problems that remain unresolved in this work.First, the low-temperature behavior, including the zero-temperature limit of the Drude weight, has not yet been calculated.In contrast, for the spin-1/2 case, the low-temperature behavior of the Drude weight has been recently addressed and calculated [68].While the low-temperature properties of the present model are more intricate than those of the spin-1/2 case, the procedure used there might be useful to our current model.
Second, the investigation of spin diffusion at the isotropic point ∆ = 1 is crucial.At this point, along with the massive region ∆ > 1, the spin Drude weight vanishes, suggesting diffusive spin transport [19,69,70].In fact, behavior similar to the divergence observed in the spin-1/2 case of the diffusion constant has also been noted in higher spin cases [19], indicating superdiffusive behavior.There is considerable interest in examining this property through more quantitative approaches, such as deriving a rigorous dynamical exponent.
Third, the Drude weight cannot be derived in our current procedure when p 0 is exactly on irrational numbers.Of course, as described before, taking the limit α → ∞, we may obtain the Drude weight when p 0 is in arbitrary vicinity of an irrational number.Indeed, the expression at infinite temperature (see (5.12)) indicates that the Drude weight converges to a certain finite value in this limit.However, it remains unclear whether this convergence value actually corresponds to the real Drude weight on p 0 when it is irrational.
Fourth, the spin Drude weight may also be evaluated using Kohn's formula [26].Namely, by imposing a twisted boundary condition on the row transfer matrix (2.20), achieved through the multiplication of the factor exp î −iϕ(n − σ z 1 )/2 ó inside the trace, where σ z 1 is the zcomponent of the spin-n/2 operator acting on V (n) 1 , the Drude weight is determined by the expectation value of the curvature of the energy spectrum with respect to the twist angle ϕ: e −βEn(0) .(6.1) In the above, E n (ϕ) denotes the nth eigenvalue of the operator obtained by taking the logarithmic derivative of the twisted row transfer matrix.Note that the integrable structure (2.21) is preserved even under the twisted boundary condition.Note that the integrable structure (2.21) is preserved even under the twisted boundary condition.Furthermore, the Hamiltonian is obtained by the same similarity transformation as in (2.26), with eigenvalues that are formally in the same form as (2.41).The effect of the twisted boundary condition is incorporated into the BAE (2.40) by multiplying the RHS with the factor e −iϕ , leading to a shift in the string center in (2.45) from λ By expanding E n (ϕ) up to second order in ϕ and explicitly evaluating (6.1), in a manner similar to Zotos's approach for the σ = 1 case, it might be possible to reproduce the formula given by (4.14).
Finally, utilizing the QTMs developed in this study, we can systematically calculate the infinite temperature Onsager matrix element associated with spin transport for the S = 1/2 chain in a zero magnetic field.This becomes feasible because the dressed scattering kernel for the S = 1/2 model at infinite temperature can be determined using the TBA equations for an S > 1/2 model.Our detailed analysis uncovers an intriguing popcorn structure of the Onsager matrix element.Further details will be reported in an upcoming publication (see [71]).
The modified TS-numbers n j defined in (A.14).
Here, v(n) is given by (2.44), and v j is the string parity defined as Using the definitions (A.15), (A.5) and (A.2), one can readily verify (A.16).Finally, we also introduce {q j } ≥1 and { q j } ≥1 as B Proof of the equivalence between (2.43) and (3.22)Let us show that the condition (2.43) is equivalent to (3.22).This is the same as showing that the set of integers n satisfying is equivalent to { n j } mα j=1 , where p 0 is a rational number given by (3.21).This can be achieved by using the method similar to that in [27].

C Derivation of the TBA equation
We briefly outline the method for deriving the TBA equations, as given in (3.34), starting from the Y -systems in (3.30).Although the procedure described in [38] can be applied heresince the Y -systems are identical to those in the spin-1/2 case (i.e., σ = 1)-the analytical properties become more complex for spins higher than 1/2 (i.e., σ > 1).As a result, the driving terms in the TBA equations manifest differently compared to the spin-1/2 case.Let As in [38], the following lemma is useful for deriving the TBA equations.
Lemma C.1.Let f j (v) for j ≥ 0 be functions satisfying where 0 ≤ v j < v 0 for j > 1.Also assume that f j (v) is analytic, non-zero, and has constant asymptotics (ANZC) in the domain v ∈ S[w j ] where v j < w j .Then f 0 (v) can be given by the following NLIE: where and the constant term is determined by the asymptotic values on both sides.
This lemma is applicable to the Y -system (3.30) after some modification to the Yfunctions defined in (3.26) and (3.27).Specifically, for m r−1 ≤ j ≤ m r − 2 (1 ≤ r ≤ α) and j = m r − 1 (r σ ≤ r ≤ α), we introduce where j σ and r σ are defined by (3.23), and the signs + and − in the exponent and before u N correspond to the sign of J, with + chosen when J > 0 (u N < 0) and − when J < 0 (u N > 0).See (3.5) for the relationship between J and u N .Additionally, for j = m r − 1 (1 ≤ r ≤ r σ − 1), we also define Again, the signs ± in the exponent and before u N in ψ r,σ (v) for j σ = m α − 1 are determined according to the specification given in (C.4).Since the identity is valid for arbitrary c ∈ C, the Y -functions in the Y -system can be replaced by the above modified ones: In the above modification, Lemma C.1 is applicable to (C.7) and the resultant NLIEs read (C.8) In the above NLIEs, the Trotter limit N → ∞ can be taken analytically and the resultant equations lead to the TBA equations (3.34).

D Free energy
Let us express the free energy, as given by (3.8), in terms of the solutions to the TBA equations (3.34).For this purpose, we first express Λ(u N , v) = T σ (v) in terms of Y -functions.Using (3.28) and (3.23), we obtain Here we have introduced the notations to save space.Referring to the definitions of the sequences in Appendix A and the periodicity (3.20), we can further modify this equation to obtain From (3.13), we observe that T n (v) possesses trivial zeros stemming from the vacuum part.Crucially, these zeros obstruct the proper application of Lemma C.1.To avoid this, we modify T n (v) to remove these zeros where m r ≤ j < m r+1 (0 ≤ r ≤ α − 1).To derive this, we have used (3.13) at the sector M = N σ/2.Thus, Comparing both sides of (D.7) as u N → 0, we conclude that the sum of the last three terms reduces to σ ℓ=1 log [ϕ(v + 2iℓ)ϕ(v − 2iℓ)] as u N → 0. The Trotter limit N → ∞ in (D.7) can be taken analytically by expanding the last three terms up to O(u N ).Using (3.23) and the properties of the sequences defined in Appendix A, these terms can be summarized as: where (D.12)

E Drude weight via GHD
Here we summarize how to derive the Drude weight from the perspective of the GHD, making sure that this paper is self-contained.The discussion here is primarily based on [7,30].
The Drude weight associated with the current densities j µ (t, x) and j ν (t, x) is described by the Kubo formula in the form of the dynamical correlation functions in the equilibrium state: where ⟨• • •⟩ c stands for the connected correlation function.(Note that D µν /β is sometimes used as an alternative definition of the Drude weight.)These current densities satisfy the continuity equation where q µ (t, x) is the charge density of the (quasi) local conserved charge Q µ (µ ∈ N).That is, Quantum integrable many-body systems with short-range interactions, including the model under consideration, possess an infinite set of conserved quantities {Q µ } ∞ µ=0 .In these systems, the maximal entropy state is, in general, characterized by the GGE, whose density matrix is given by where β µ represents the generalized inverse temperature.Let us denote as the average of any operator o over the GGE.In particular, let q µ and j µ be the average charge and current, respectively: Then, the relations ∂ β µ q ν = ∂ β ν q µ and ∂ β µ j ν = ∂ β ν j µ hold [7], which ensure the existence of the free energy density and the free energy flux g (whose explicit form is defined later) such that Now, we introduce the equal-time charge-charge and current-charge correlation functions: q µ C µν (q ν , o). (E.14) The matrix C µν denotes the inverse of C µν .Thus, as referenced in [30], the Drude weight can be rewritten as β −1 D µν (β) = γ,δ (j µ , q γ )C γδ (q δ , j ν ).(E.15) Using the definition (E.9) and the relation (E.12), we find that D µν (β) is given by the following matrix forms: Let us now apply the argument discussed above to a generic fermionic quantum integrable system (which of course includes the present model), whose thermodynamic quantities can be described by the TBA.In such a system, the local conserved charges q µ are decomposed in terms of the quasi-particles of type j (corresponding to the n j -string in the case of our present model): q µ = j dv ρ j (v)h µ;j (v) = j dv ρ tot j (v)ϑ j (v)h µ;j (v).(E.17) Here, ρ j (v) is the distribution density of the quasi-particles of type j, and ρ tot (v) := ρ j (v) + ρ h j (v) denotes the density of states, where ρ h j is the density of holes.The quantity h µ;j (v) is referred to as a bare charge.In particular, we set to be the bare energy, and accordingly, The quantity represents the fermionic distribution function.From the BAE (see (2.40) for instance), we obtain ρ tot j (v) as the solution to the integral equation: where ξ j is a sign factor, and K jk (v) is a kernel assumed to be symmetric K jk (v) = K kj (v).The quantity p j (v) is the bare momentum, related to the bare energy (E.18) by p ′ j (v) = ϵε j (v) = ϵh 1,j (v), (E.22) where ϵ ∈ R is a constant depending the definition of the overall scaling factor of the Hamiltonian.For instance, ϵ = − θ J sin θ (E.23) for our present model.In the framework of the TBA, the free energy density f (E.7) is written as Correspondingly, the free energy flux g (see (E.8)) is given by [5] as The function η j (v) defined by The above dressed operation can also similarly apply for any function h j (v) by just replacing h µ,j (v) with h j (v).Some crucial quantities related to the Drude weight are expressed by the dressed functions.For example, the density of state ρ tot j (v) (E.21) is expressed as the dressed energy: This can be easily followed from (E.21) in conjunction with (E.22) and (E.29).In particular, the effective velocity v eff j (v) of the quasi-particle of type j is given by In addition to (E.17), we derive q µ using (E.8) and (E.24): q µ = 1 2π j dv p ′ j (v)ξ j ϑ j (v)h dr µ;j (v).(E.32) Noticing the relation in (E.22) and applying the dressed function formalism to (E.29), we find that (E.32) aligns with (E.17).Similarly, from (E.8) and (E.25), we obtain j µ = 1 2π j dv ε ′ j (v)ξ j ϑ j (v)h dr µ;j (v) = 1 2π j dv ∂ v (ε dr j )(v)ξ j ϑ j (v)h µ;j (v) = j dv ρ j (v)v eff j (v)h µ,j (v), (E.33) where the dressed function formalism is applied to derive the second equality.Eq. (E.30) has been used to derive the third equality which provides a physically reasonable picture.The relation (E.33) was originally proposed by [5,6], and has recently been confirmed by the Bethe ansatz technique [16,75,76].According to the definitions given in (E.9), we differentiate (E.32) and (E.The relation between B and C as given by (E.12) (more explicitly B µν = γ A γ µ C γν in terms of the matrix elements) leads to A ν µ h dr ν;j (v) = v eff j h dr µ;j (v).(E.36) This indicates that the dressed charges serve as the eigenfunctions of the flux Jacobian, and their corresponding eigenvalues are the effective velocities.Then, using (E.16) with (E.36), we finally arrive at ó 2 h dr µ;j (v)h dr ν;j (v).(E.37) All the quantities in the integrand are calculated using the generalized TBA equations (E.27).
After replacing them with those obtained via the standard TBA, we also evaluate the original Drude weight D µν (β) (E.1) by (E.37).
Therefore, the O(β) contributions from ϵ k l to T (u N , v) cancel each other out.Conversely, the contributions from δ k l are retained in the form of the total sum δ k := N/2 l=1 δ k l .The quantity δ k is obtained as follows.First, substituting u = u N and (3.15) into (3.14), and using the fact that z k l := ϵ k l + iδ k l and u N are both O(β), we obtain