Thermodynamic Bethe Ansatz and Generalised Hydrodynamics in the sine–Gordon model

We set up a hydrodynamic description of the non-equilibrium dynamics of sine–Gordon quantum field theory for generic coupling. It is built upon an explicit form of the Bethe Ansatz description of general thermodynamic states, with the structure of the resulting coupled integral equations encoded in terms of graphical diagrams. The resulting framework is applied to derive results for the Drude weights of charge and energy. Quantities associated with the charge universally have fractal dependence on the coupling, which is notably absent from those associated with energy, a feature explained by the different roles played by reflective scattering in transporting these quantities. We then present far-from-equilibrium results, including explicit time evolution starting from bipartite initial states and dynamical correlation functions. Our framework can be applied to explore numerous other aspects of non-equilibrium dynamics, which opens the way to a wide array of theoretical studies and potential novel experimental predictions.

Various methods can address non-equilibrium time evolution in sine-Gordon quantum field theory.The full quantum field dynamics can be addressed using truncated Hamiltonian approaches.However, these are limited in time by finite volume and to small quenches by the presence of the cutoff [27,28].Another approach to quantum dynamics uses the exact form factors of local operators [29,30] to construct a spectral expansion valid for small quenches.However, its extension to the attractive regime [31] faces severe difficulties [32], which remain unresolved.
Semiclassical field theory approximations, such as the truncated Wigner approximation [33] or a self-consistent Hartree-Fock method [34,35] are limited in scope to the deep semiclassical regime and contain approximations which are hard to control [27,28,36].A different semiclassical approximation [37][38][39][40] can be constructed using the quasi-particle picture of time evolution.Although the quasi-particle picture is more generally valid for time evolution induced by quantum quenches [41,42], the semiclassical approximation can only take into account a very simplified version of the scattering processes, although it can be supplemented by some quantum corrections [39].
Description of the dynamics in inhomogeneous situations is much more limited, with only some approaches having extensions to this case, such as the semiclassical quasiparticle approach [40] and truncated Hamiltonian methods [43].
Recently, a novel approach has been developed to describe the large-scale non-equilibrium behaviour of integrable systems called generalised hydrodynamics (GHD) [44,45].It relies on the separation between the scales of spatial variation and microscopic quantum dynamics scales, which occurs in many physical situations.This leads to the assumption of local thermodynamic equilibrium, a common basis for hydrodynamic approaches in general.The specific feature of integrable systems is an infinite number of local conserved quantities, leading to an infinite set of continuity equations.To obtain a closed system of equations, it is necessary to describe the locally homogeneous thermodynamic state, which for integrable systems is accomplished by the thermodynamic Bethe Ansatz (TBA) [46][47][48].However, for the sine-Gordon model, the main issue preventing progress in modelling the non-equilibrium dynamics has been the absence of an explicit description of thermodynamic states for general couplings.So far, it has only been formulated explicitly for special values of the coupling [49][50][51], although a corresponding set of functional relations (the so-called Y system) was conjectured for the general case in [52].We note that the thermodynamic description of the classical sine-Gordon model has only been constructed recently [53].
This work explicitly derives the thermodynamic Bethe Ansatz system necessary for formulating generalised hydrodynamics for sine-Gordon theory at general coupling values.We then set up the sine-Gordon GHD and consider its predictions for transport in the theory.In particular, we compute Drude weights for energy transport, comparing them to previously obtained charge Drude weights [54], and also examine the asymptotic state obtained after a bipartition protocol.
The outline of the paper is as follows.In Section 2 we present the sine-Gordon TBA system for thermal states in a partially decoupled form following the example of the XXZ spin chain [48].We obtain the ingredients required for the GHD in Section 3, which include equations for quasi-particle densities, the so-called dressing relations and effective velocity.In Section 4 we study the charge and energy Drude weights.Section 5 presents results obtained from the GHD, such as the asymptotic state resulting after time evolution in a bipartition protocol and dynamical two-point correlation functions.We present our conclusions and outlook in Section 6.The detailed derivations of the coupled and partially decoupled TBA systems are presented in the Appendices.

Thermodynamic Bethe Ansatz of the sine-Gordon model
This section presents the TBA equations for generic couplings in fully coupled and partially decoupled forms and tests them against the Destri-de Vega NLIE.

The sine-Gordon model and its factorised scattering theory
The dynamics of the sine-Gordon field theory is governed by the Hamiltonian (we set ħ h= c =1) where φ is a real scalar field, and normal ordering is defined relative to the modes of the free massless boson obtained in the limit λ = 0.The coupling λ is dimensionful and defines the scale of the theory, with the eventual strength of interaction determined by the dimensionless coupling β, which takes values 0 < β 2 ≤ 8π for which the cosine interaction is relevant.
To describe the spectrum and the exact scattering theory, it is convenient to introduce the renormalised coupling The spectrum contains a doublet of topologically charged excitations consisting of a kink and antikink of mass m S , which is related to the parameter λ as [55] λ = . (3) In the repulsive regime ξ > 1, the kink doublet is the only particle excitation in the spectrum, while in the attractive regime ξ < 1, there are ⌊1/ξ⌋ species of neutral kink-antikink bound states also known as breathers with masses except for the points where ξ −1 is integer, where the number of breathers is 1/ξ − 1.The momenta and energies of these particles are given in terms of the rapidity θ as p(θ ) = m a sinh(θ ), E(θ ) = m a cosh(θ ).
The theory is integrable at the classical and quantum levels, implying that general multiparticle scattering processes can be factorised in terms of 2 → 2 particle processes.The scattering of the kinks is described by the following two-particle amplitudes [15]: where +/− stands for kinks/antikinks, with θ denoting the difference of their rapidities.Note that for integer values of 1/ξ, the kink-antikink reflection amplitude S R vanishes, corresponding to purely transmissive scattering; these points in the parameter space are called reflectionless.
The breather-soliton and soliton-soliton scattering amplitudes can be specified in terms of the following elementary blocks: Since these amplitudes are pure phases, their logarithms define phase shifts, which can be fixed unambiguously by a suitable choice of the branch of the logarithm.We adopt a convention for which the phase shift δ a corresponding to an elementary block is defined as with δ a (0) = 0 and δ a (±∞) = π.The derivative of the phase shift is Throughout this work, we adopt the following conventions for Fourier transformation and the convolution which means that the Fourier transform of ϕ a can be written as With these preliminaries, the scattering amplitudes between a kink/antikink and a breather are while the amplitudes for the scattering of two breathers are All entries except the first and the last are repeated, indicating that the corresponding block occurs twice, i.e.P nm is treated as a set with multiplicities (multiset).

The sine-Gordon TBA system
Here we discuss the sine-Gordon TBA system apart from reflectionless points. 1 For the details of the derivation, see Appendix A.
Let us write the coupling as with n B specifying the number of breathers and α ≥ 1, with reflectionless points corresponding to setting α = 1.Considering the field theory in a finite volume R, we can describe a generic eigenstate of the system by specifying the following ingredients: • Rapidities of breathers θ , where k = 1, . . ., n B runs through the breather species and for a fixed k, j = 1, . . ., N B k where N B k is the number of breathers of species k present; • solitonic rapidities θ ( j) S , j = 1, . . ., N S , where N S is the total number of kinks/antikinks present; and • so-called magnonic rapidities that are variables specifying the internal wave function in the 2 N S -dimensional space of topological charges.The number of species (breathers, one or two solitons, and magnons altogether) appearing in the Bethe Ansatz equations as a function of the coupling strength shows an intricate structure.We adopt the convention of denoting couplings corresponding to zero (reflectionless), one, two, and three magnonic levels with green circles, blue squares, red diamonds and brown stars, respectively.
Up to exponentially small corrections in the volume, the energy of the state (relative to the vacuum) can be computed as where the rapidities satisfy the Bethe equations (A.1).The magnons can be brought into one-to-one correspondence with the Bethe Ansatz of the gapless XXZ spin chain of length N S , which allows us to borrow the string hypothesis for the magnons from the XXZ spin chain: the thermodynamic limit is assumed to be dominated by string configurations of elementary magnons which we call magnonic strings or simply magnons.The string configurations can be classified writing ξ as a (unique) simple continued fraction where the ν i , i = 1, . . ., l are positive integers.We call the number l of integers ν i appearing in the above representation the number of levels, which is finite when ξ is rational and infinite when it is irrational.The number of magnonic species is given by the sum of the integers, and they can be indexed with a species label M i with i = 1, . . ., n M .As shown in Fig. 1 for several rational couplings, the number of species has a very complicated dependence on the couplings.Magnonic strings (or magnons for short) consist of elementary magnons with the same real and equidistant imaginary parts.The number of elementary magnons that make up a string is called the string's length, which we denote by ℓ.In the thermodynamic limit, the Bethe Ansatz is formulated in terms of the densities of Bethe Ansatz roots (filled states), denoted by ρ i (θ ).

The difference ρ
is the density of unoccupied rapidities called holes.
We find where the integration kernels Φ are given explicitly in Appendix A in Eqs.(A.19,A.22).Observe that the system (16) has the overall form where m i = 0 for magnonic degrees of freedom.Note that the breathers and the magnons are coupled only through the soliton, not directly.Following the usual procedure [46], the TBA equations for the thermal equilibrium state follow by minimising the free energy density with respect to the root densities ρ i subject to the conditions (17).Here T is the temperature, s is the Yang-Yang entropy density [46], and µ is the chemical potential coupled to the topological charge, while the q i are the topological charges carried by the various excitations: Introducing the pseudo-energy functions the resulting TBA system is which can be written in the concise form where the source terms are w i = m i cosh θ /T − µq i /T .The free energy density f of the equilibrium state can be computed as

Partial decoupling of the TBA system
As noted in the previous section, the magnonic part of the sine-Gordon TBA is essentially identical to the TBA system of the XXZ spin chain, for which a partial decoupling of the system of equations can be achieved [56], which is also the case for quantum field theories with diagonal scattering [49].Partial decoupling means that a system like ( 21) can be recast in the form where every pseudo-energy is only coupled to a few others, and the structure of the TBA can be represented with a simple graph.Moreover, the kernels in this form are also much simpler, and the partial decoupling makes the system's numerical solution much more computationally efficient.
Albeit the decoupling problem is solved for the XXZ spin chain, the result cannot be directly transferred to the system (21) due to the presence of the massive nodes, and the structure must be analysed carefully.We adopt the following strategy.First, we decouple the system with a single magnonic level, i.e. when the coupling has a continued fraction expansion It turns out that the major issues can be fixed by considering this case.We then perform decoupling with two magnonic levels, i.e. when We find that the resulting structure is stabilised with the equations involving only magnons coinciding with the appropriate decoupled equations for the XXZ spin chain using the mapping defined by Eqs.(A.4),(A.5).It is then clear that higher magnonic levels can also be directly obtained from the XXZ case, and we check this by considering systems with three and four magnonic levels.The validity of these systems can be checked numerically in a very stringent way by comparing the free energy at a finite temperature and vanishing chemical potential to that resulting from the Destri-de Vega (DdV) complex nonlinear integral equation [57].In addition, we also carry out various self-consistency checks.
After decoupling, the TBA system can be concisely written as where L j = log 1 + e −ε i and the kernel K i j is a sparse matrix encoding the coupling of each species to the others, σ (1,2) j are numerical coefficients taking values 0 and 1, and the source Table 1: Building blocks of diagrams encoding the sine-Gordon TBA systems at different couplings.
terms w i are modified by the decoupling procedure in comparison to the source terms w i appearing in the fully coupled TBA equations (22).This modification does not affect the contributions related to energy and momentum; however, for all other charges, their oneparticle value h i (θ ) is modified to a new form h i (θ ).In the cases considered here, this only affects the topological charge.However, when extending the TBA system by higher conserved charges to construct a generalised Gibbs ensemble [58], all the corresponding one-particle eigenvalues are modified in the decoupled equations.
For later reference, we also introduce a notation for the kernels that appear in this section specified by their Fourier transforms with the p i defined in Eq. (A.10).
To depict the kernel matrix K i j in Eq. ( 26) in a graphical way, we borrow a formalism introduced in Ref. [59] which treated the TBA of the boundary sine-Gordon model.We extended this formalism to allow for special cases when some of the integers in the continued fraction expansion (15) are small, i.e. satisfy ν i ≤ 2. The diagrams can be built using six types of blocks corresponding to various kernels as summarised in Table 1.
The decoupling procedure up to two levels is detailed in Appendix B, while the resulting TBA system at general coupling is presented in the next section.

The TBA system at general coupling
Once we have the system for one and two magnonic levels, all possibilities for sewing together the massive and the magnonic nodes through the soliton are covered.Deeper magnonic levels are exactly identical to the XXZ Bethe Ansatz, so instead of explicitly performing the decoupling, the results can be borrowed from that case [48].Note that the system's structure is somewhat altered at reflectionless points, for which the result is given in Eq. (B.3).
Introducing the notation Table 2: Source terms and coefficients appearing in the TBA system Eq.( 26) and the dressing equation ( 47) in the generic case (except the reflectionless points).See Eqs.(A.10,A.12)for the definition of y i and r i .This table extends the results of [54], which were only derived for up to two magnonic levels.
the resulting system takes the form where Θ is the Heaviside function and δ is the Kronecker delta, with the parameters appearing in the above system defined in Eqs.(A.10,A.11,A.12)and summarised in Table 2. To ease the notations further, we introduced a species label M 0 whose occurrences should be replaced by the soliton label S on the RHS, while equations with M 0 on the LHS must be omitted.
Appendix B shows examples of the graphical representation of Eq.( 29).Figs. 2 and 3 depict more intricate graphs with multiple magnonic levels with a different number of magnons on each level to show examples with various structures.Fig. 3 also shows how a graph is altered as the number of magnons ν 2 at level 2 changes.
So far, we assumed that the continued fraction (15) has a finite number of levels corresponding to a rational value of the coupling parameter ξ.We note that irrational couplings correspond to an infinite continued fraction (15); truncating the continued fraction at progressively deeper levels leads to rational approximants of the coupling converging to the eventual irrational value.Similarly to the case of the XXZ spin chain, the relevant physical quantities are expected to be obtained as a limit of a sequence constructed from these rational approximants.This is trivially valid for quantities (such as the free energy) that are a smooth function of the coupling.For quantities with a fractal dependence of the couplings, such as charge Drude weights, numerical evidence shows that the discontinuities decrease with the number of magnonic levels and are fully consistent with convergence in the limit of the infinite continued fraction.

Ultraviolet limit and central charge
In the high temperature (ultraviolet) limit, the cosine perturbation in the sine-Gordon Hamiltonian (1) can be neglected, and the theory becomes that of a massless boson corresponding to the first two terms.We can test our TBA equations by checking that they are consistent with this behaviour.
In the limit of large temperature, T ≫ m S , m B k , the source terms in the TBA equations (29) of the massive particles of mass m a can be neglected in the rapidity domain − log(2T /m a ) ≲ θ ≲ log(2T /m a ).Consequently, all the source terms are constant (in a thermal state), so the pseudo-energies also take constant values ε a , and the TBA equations reduce to a set of algebraic equations for the ȳa = exp( εa ) "plateau" values.
For example, for two magnonic levels, using that the integral of all kernels Φ a are equal to π (equivalently, Φ(0) = 1/2), we obtain in the generic case , demonstrating how the TBA system collapses at level 2 as the corresponding number of magnons is gradually reduced. where For µ = 0 the solution for two magnonic levels is For a single magnonic level, the same expressions hold up to magnon number ν 1 − 2, and Interestingly, these expressions also hold in the special cases (e.g.ν 1 = 1 or ν 2 = 2) and in the repulsive regime with n B = 0.In the reflectionless case, there are no magnons, but there are two soliton nodes instead of one (the soliton and the antisoliton).The breather pseudoenergies are the same as in the generic case, and the soliton and antisoliton plateau values are y S = y S = n B + 1.
In the complementary domain |θ | ≫ log 2T /m a , the pseudo-energy functions for the massive excitations grow rapidly as ∼ m a cosh θ and their filling functions decay to zero faster than exponential.As e ε S (θ ) becomes very large, it can be dropped in the equation of ε M 1 (θ ).As a result, the magnonic TBA equations decouple completely from the massive equations.Since the source terms are constant in a thermal state, the magnonic pseudo-energies become constants ε M k , and y M k = exp( ε M k ) obey algebraic equations again.These are the same as the magnonic equations in (30) except for the first magnon, which reads The solution for two magnon levels and µ = 0 reads For a single magnonic level the first line holds and The root densities ρ a (θ ) are supported around θ = log 2T /M a and exponentially small elsewhere, so the equations split to equations describing independent left and right moving modes.In these rapidity regions, the effective velocities agree with the speed of light ±1 for all the excitations.
The free energy density can then be expressed using the standard procedure [47, 60, 61] based on Roger's dilogarithm function Taking into account the nonzero constant pseudo-energies of the magnons as |θ | → ∞ as well as their associated signs η M k , the free energy density in a thermal state with µ = 0 is where are the constant filling fractions.Substituting the solutions of Eqs. ( 34) and ( 31), we checked numerically in various cases that which is the exact result for a free massless boson (a conformal field theory with central charge c = 1).This provides another consistency check of the validity of our TBA equations.We note that, from a mathematical viewpoint, Eq. ( 37) constitutes nontrivial dilogarithm identities [52,62].
Figure 4: Comparison of the free energy density calculated from DdV (Eq.( 41), black curve) and from TBA (Eq.( 23), symbols) for two temperatures and several different values of the coupling, shown with coloured symbols as specified in Fig. 1.The relative difference between the result of the two methods is less than 10 −5 in all cases we considered.

Comparison with the NLIE
An independent verification for our TBA system can be obtained by setting the chemical potential to zero and comparing the resulting free energies to those computed from the Destri-de Vega complex nonlinear integral equation (NLIE) [57] The above equation can be solved iteratively for the function Z(θ ), from which the free energy density f is obtained using the formula Comparisons of the free energy density calculated from the DdV and the TBA for different values of the coupling and temperature are shown in Fig. 4 and Tables 3, 4. The excellent agreement provides a nontrivial and stringent verification of the TBA system (29).

Dressing and further tests of the TBA system
This section presents the dressing equations and the partially decoupled form of the density equations.We also perform further consistency checks of the full TBA system.

Dressing and partially decoupled equations for the densities
The presence of finite quasi-particle density in thermodynamic states leads to a dressing of all one-particle quantities, such as momentum, energy, and charges.To derive the dressing equations, we follow [63] and write the source terms in Eq. ( 26) as where h i (θ ) are the charge eigenvalues modified by the decoupling procedure.The oneparticle energies and momenta e i = m i cosh θ and p i = m i sinh θ are unchanged, while q i are the topological charge values resulting after decoupling which must be distinguished from the bare charges q i .The "partially decoupled charges" are given by for massive particles, and intermediate magnons, − y l , for the last two magnons, (43) as can be seen from Table 2.Note that in the reflectionless case, the charge assignment is In addition, β (e) = 1/T , β (p) = λ/T , and β (q) = µ/T are the thermodynamic conjugate variables associated with energy, momentum and topological charge.We note that this idea can be extended to construct generalised Gibbs ensembles from TBA by including the higher conserved charges associated with integrability [58].
Starting from the free energy (23), the expectation value of a charge conjugate to the generalised temperature variables β (h) where h dr i is the dressed charge, and we introduced the filling fractions: Using Eqs.(26,42), the dressed charges satisfy the dressing equation In particular, the total density of states corresponds to dressing the derivative of the momentum (which is equivalent to the energy), i.e.
These equations are nothing but the partially decoupled versions of Eqs.(A.24), which formulate the Bethe Ansatz (A.16) for thermodynamic states in terms of the densities.We stress that these equations hold for all thermodynamic states, i.e. for states in which the quasi-particles can be described in terms of density functions.The density equations provide relations between the total densities of states ρ tot i and the filling fractions ϑ j (or, equivalently, the root densities ρ i ).Still, they do not determine the state by themselves.In thermal equilibrium, the missing information is provided by the pseudo-energy functions which solve the TBA equations Eqs.(A.28) or, equivalently, their partially decoupled form Eqs. (29).In an inhomogeneous situation, the missing information is provided by the time evolution determined by the GHD equations (60) introduced in Section 5.
Comparing Eq. ( 48) to (47) demonstrates that the density equations can be obtained by taking the derivatives of the TBA equations [48], which can be compared to the result of explicitly decoupling the density equations, which is eventually used in Subsection 3.2 as a cross-check for our calculations.
The dressed values of the (rapidity derivatives of) energy and momentum can be used to define the effective velocity which enters the equations of Generalised Hydrodynamics (60) and can be interpreted as the velocity of individual quasi-particle excitations in the finite density medium constituted by the other quasi-particles present in the thermodynamic state.

Further tests of the TBA system
In the previous subsection, the description of thermodynamic states was extended by introducing dressing, an example of which is Eq. ( 48) for the densities.In Figs. 5 and 6, we show examples of what the fillings, the root densities, the effective velocities and the dressed topological charges look like as functions of the rapidity.We found that in the attractive case, a shorter θ -grid is enough to resolve the non-trivial structure of the above quantities, while in the repulsive case, this grid usually needs to be longer.We tested the self-consistency of the TBA system and our numerical implementation by: (i) Comparing free energy values calculated from the coupled system (A.29) and the partially decoupled system (26).
(ii) Comparing free energy values calculated from the pseudo-energies ( 23) and the densities (18).
(iii) Comparing numerical derivatives of the pseudo-energy functions as in the derivation in Sec.3.1, to the direct calculation of (47).This also provides a way to check the signs η i , as they do not appear in the pseudo-energy equations ( 26) but do appear in the dressing equations (47).
(iv) Testing the symmetry of the free energy under the µ → −µ transformation, which must hold on physical grounds but is not manifest in the TBA system (26).
(v) Calculating charge and current expectation values in multiple ways, which must give Note that here h i is the bare charge and h dr i is its dressed counterpart, but as mentioned in connection with Eq. ( 42) above, in the partially decoupled form of the TBA equations ( 26) their contribution to the source term is modified to h i .In each case, we have found good agreement, confirming the validity of the TBA system and its self-consistency.

Transport and Drude weights
As our first application of sine-Gordon thermodynamics, we consider Drude weights describing ballistic transport of charge and energy.

Drude weights from TBA
The Drude weight is defined as the integral of the connected correlation function: where j h is the current of the charge, which we specify by giving its one-particle eigenvalue h i (θ ) that depends on the species i and rapidity θ of the excitation.In the TBA framework, the Drude weight can be computed as [64-66] where the sum runs over all particle species.However, for topologically neutral states (µ = 0), the dressed topological charge of all species turns out to be zero except for the last two magnons with opposite and constant dressed charges.
The energy and charge Drude weight results are shown in Fig. 7.The charge Drude weights show the characteristic fractal pattern already established in our previous work [54], while the energy Drude weights depend on the coupling in a fully continuous manner.This is due to a fundamental difference in energy and charge transport.In the scattering of field-theoretic excitations (kinks and breathers), the energy (parameterised by rapidity) is always transmitted.However, in the kink-antikink scattering described by the amplitudes (5), the topological charge can also be subject to reflection, altering its transport.The non-diagonal structure of the kink-antikink scattering is reflected in the magnonic quasi-particles of the Bethe Ansatz, which have a very intricate structure that follows the continued fraction representation (15) of the coupling.The observed fractal structure is similar to the case of the XXZ chain in its gapless regime [67][68][69][70][71], a phenomenon also known as "popcorn" Drude weights [72].The sine-Gordon model is the second example for which the fractal nature of the spin Drude weight was established; more recently, it was established for higher spin variants of XXZ spin chain as well [73].
It is argued in the recent work [72] that the fractal structure of the Drude weight results from the U q (sl(2)) symmetry of the model.This argument is strongly supported by finding fractal Drude weights in the sine-Gordon model since it has the same quantum symmetry [74,75], which is also true for the higher spin variants of XXZ spin chain.Note that while the magnonic part of the Bethe Ansatz is essentially the same as that of the gapless XXZ spin chain, the presence of the massive breather and solitonic excitations of the Bethe Ansatz makes the sine-Gordon TBA system essentially different from that of the spin chain.Therefore, the finding of a fractal structure in the sine-Gordon model, while not unexpected, is still a nontrivial confirmation of the arguments advanced in [72].
Note that in a naive semiclassical picture, the presence of reflections results in a random walk for the topological charge, which would seem to imply a diffusive behaviour.Therefore, the nonzero charge Drude weight requires a specific explanation.In the framework of Mazur's inequality [76,77], a nonzero Drude weight for the topological charge strongly suggests the existence of yet unknown conserved quantities that are odd under charge conjugation, in parallel with those found in the XXZ spin chain [68].

Low-temperature limit at the reflectionless points
In the attractive regime and at reflectionless couplings ξ = 1/(n B − 1), the TBA simplifies to a form involving only breathers and a kink/antikink pair given in (B.3).For low temperatures, the convolution terms in the TBA equations can be dropped, and all effects of interactions are exponentially suppressed due to the mass terms.As a result, the kinks/antikinks can be described as non-interacting fermions with energy M cosh θ and with an effective velocity equal to their bare velocity tanh θ .The filling factors are given by simple Fermi-Dirac distributions    Even though the number of particles has an intricate dependence on the coupling (cf.Fig. 1), the TBA system captures the continuous coupling dependence of the energy Drude weight.In contrast, the charge Drude weight shows a fractal/popcorn dependence on the coupling strength.Besides the numerically calculated values, the bottom left panel shows Eq. ( 54) with a red dashed line, showing excellent agreement at the reflectionless points.In the rightmost figures, we omit individual points for clarity and only display the data using continuous black lines, as the fractal structure is not resolved numerically at these high temperatures.The continuous behaviour of the charge Drude weight (bottom right panel) in the limit T → ∞ is justified by the agreement with the result (55) of the analytical calculations in Sec.4.3, shown as a red dashed line.The analytical high-temperature limit of the energy Drude weight, Eq. ( 57), plotted as a red dashed line in the top right panel, also shows good agreement with the numerical results.
Due to the absence of interactions, the dressing ( 47) is trivial since all kernels K i j vanish.Additionally, all relevant signs are trivial (η i = +1); therefore, all quantities are equal to their dressed counterparts.Since only the kinks carry topological charge, the Drude weight (52) simplifies to the explicit expression with the factor 2 accounting for the presence of kinks and antikinks, which carry topological charge h(θ ) = h dr (θ ) = ±1.The result ( 54) is independent of the coupling and agrees fully with the numerical results obtained from the full TBA as shown in Fig. 7.

High-temperature limit
In the high-temperature limit, the sine-Gordon interaction can be neglected, and the dynamics can be described by the Hamiltonian of a massless free boson, which makes possible the explicit evaluation of the charge Drude weight with the result We refer to the Supplemental Material of [54] for details.This result agrees perfectly with the numerically obtained data, as shown in Fig. 7.Note that the high-temperature limit of the Drude weight is a continuous function of the coupling parameter ξ with the fractal structure suppressed, except when β 2 is close to 8π.In fact, the numerical computations described in the previous subsection show that the Drude weight goes to zero at all finite T in the Kosterlitz-Thouless limit β 2 /8π → 1; therefore Eq. ( 55) implies that the limits β 2 /8π → 1 and T → ∞ do not commute.
The high-temperature limit of the energy Drude weight can be obtained simply from known results in non-equilibrium conformal field theory.In a bipartitioned system with temperatures T 1 /T 2 in the left/right halves [78], respectively, the current flowing between the two halves is where c is the central charge, which in our case is 1.Using the formula (65) below results in which agrees well with the numerics as shown in the top right panel of Fig. 7.In addition, the energy Drude weight also shows the non-commutativity of the limits β 2 /8π → 1 and T → ∞.

Charge Drude weight at finite chemical potential
To further analyse the fractal structure of the charge Drude weight, we also calculated it for finite chemical potential.The results, shown in Fig. 8, reveal that the fractal structure persists away from zero chemical potential but is gradually suppressed with increasing µ.

Generalised Hydrodynamics in the sine-Gordon model
Generalised Hydrodynamics is a framework that describes the dynamics of integrable systems on the hydrodynamic (Euler) scale.It is based on the transport of the infinitely many conserved quantities captured by the continuity equation The expectation values of the charge and current density can be expressed in terms of the root densities as The expression for the current densities was originally conjectured in [44,45] and proven later in [79].Exploiting the completeness of the charges, one arrives at the GHD equation [44,45]     for the densities ρ i (t, x, θ ) of quasi-particles of species i and rapidity θ that are space and time-dependent on the Euler scale.The effective velocity (49) carries an implicit dependence on t and x via the densities {ρ j } used to dress the derivatives of energy and momentum.These equations are supplemented by the dressing equations (47) and the density equations ( 48) necessary to reconstruct the filling fractions needed for the dressing from the quasi-particle (root) densities ρ i .

Bipartition protocol
Arguably the most frequently implemented protocol to study non-equilibrium dynamics in inhomogeneous states is the bipartition protocol, where the two halves of an infinite system are prepared in different equilibrium states and the system is described by the source terms At t = 0 the system is let to evolve freely, and it is shown e.g. in [44,45] that in the limit x, t → ∞, the state of the system is described by the filling function where ζ = x/t and ϑ i,L and ϑ i,R are the filling functions of the left and the right initial states.Note that Eq. ( 63) is an implicit equation for ϑ i , as the effective velocities on the right-hand side depend on ϑ i .Despite this, we found that the usual recursive scheme [63,80] also converges to the solution in the sine-Gordon model.Using Eq. ( 50) together with the fillings at each ray, one obtains the asymptotic h(ζ) and j h (ζ) profiles in the bipartite system, which are the limits along "rays" of fixed ζ = x/t: Examples of such energy and topological charge profiles are shown in Fig. 9.Note the cusps [81] in the charge and charge-current profiles in the repulsive regime, which are due to the maximum magnonic velocities being considerably smaller than the speed of light (corresponding to the speed of sound in a condensed matter context and equal to 1 in our units), cf.Fig. 6.In contrast, there are no such cusps in the attractive regime nor in the energy profiles in both regimes because the maximum soliton and magnon velocities are equal to the speed of light in these cases, cf.Fig. 5.
One can compare the position of these cusps to the maximum values of the magnonic effective velocities in the "reference" thermal states corresponding to the post-quench temperature T = 2m S and chemical potential µ = 0 shown in Figs. 5, 6.Note that while the maximal effective velocities qualitatively agree with the locations of the cusps, the eventual numerical values of ζ where the cusps occur in the non-equilibrium evolution differs from those one would guess from the equilibrium effective velocity profiles; for example, they are not symmetric under ζ → −ζ.The reason is that in the non-equilibrium evolution induced by the bipartition protocol, there is a different asymptotic state at each ray ζ, which also differs from the reference thermal equilibrium state.
The Drude weight of any conserved charge can also be evaluated from an infinitesimally unbalanced bipartition protocol using the linear response formula [69,82] which gives identical results to the method used in Sec. 4 [54].

Dynamical correlators
The GHD framework provides access to dynamical correlation functions on the Euler scale in a large distance/time limit.Here, we consider the simplest kind, the connected correlators of conserved densities.In a homogeneous equilibrium state, these are given by [64] where ζ = x/t and θ * i (ζ) is the set of rapidities where the effective velocity takes the value of ζ, i.e. the solution of the equation Note that the sharp peaks at higher temperatures (top two rows) correspond to the maximal value of the effective velocity, c.f. Figs.5,6.In contrast, at low temperatures, the peaks appear only when this maximum value is small (bottom right panel).Otherwise, at low temperatures, high rapidity states (corresponding to the maximum value of the effective velocity) are unoccupied; therefore, the correlator takes a bell shape.This formula expresses that correlations are built by quasiparticles that travel ballistically at velocities v eff i and contribute with their dressed charges.The correlators of the energy density and the topological charge density are shown in Fig. 10 for the attractive (ξ = 2/7) and the repulsive (ξ = 7/2) regime.For high temperatures (top two rows), they are seen to be strongly peaked at the boundaries of their support in ζ, which for the energy density corresponds to the light cone ζ = ±1.The reason for the peaks is that high rapidity states are also occupied at high temperatures, which carry correlations with the maximum of the corresponding effective velocity.For small temperatures (bottom two rows), the correlators often take a bell shape because the density of occupied states is concentrated at small rapidities, and therefore the particles cannot sample the maximal value of the effective velocity.If the maximum velocity is also small at low temperatures, the correlators take a peaked shape again, as the bottom right example shows in Fig. 10.Interestingly, the topological charge density correlator has a smaller support in the repulsive regime, indicating a difference between charge and energy transport, which we address in a separate study [83].
These correlators can also be used to reconstruct the Drude weights as [64] however, evaluating the integral accurately requires a sufficiently large number of grid points in ζ due to the presence of cusps and is very costly in numerical terms.For this technical reason, a suitably accurate calculation only proved possible in the repulsive regime, where we found complete agreement with the Drude weights computed from Eq. ( 52).

Conclusions and Outlook
In this work, we described the thermodynamics and hydrodynamics of the sine-Gordon quantum field theory at generic couplings in detail.Thermodynamic states can be specified with quasi-particle densities in rapidity space, which satisfy linear integral equations derived from the Bethe Ansatz.These equations introduce a relation between the total and occupied (root) densities of states, allowing us to determine one set of densities in terms of the other.For equilibrium states such as grand canonical (or arbitrary generalised Gibbs) ensembles, a system of nonlinear integral equations known as Thermodynamic Bethe Ansatz (TBA) supplies the missing information to determine the state completely.For states inhomogeneous at the hydrodynamic Euler scale, the Generalised Hydrodynamics (GHD) equations determine the evolution of the root densities from an arbitrary initial condition given in terms of the root densities.GHD requires the determination of the effective velocity of quasi-particles as an input from the Bethe Ansatz, which is given in terms of the dressing equations.After deriving the full thermodynamic description, including the density equations, the TBA, and the dressing equations, we verified their structure extensively, cross-checking them against each other and comparing them to the Destri-de Vega nonlinear integral equation.We presented the thermodynamic equations in both their fully coupled and partially decoupled form, and the derivation of the latter was greatly simplified by incorporating results from the XXZ spin chain.The decoupled equations have several advantages, such as (i) the much simpler structure of kernels and (ii) the sparse form of coupling between densities, which greatly reduce the computational cost associated with their solution.This is very important in the GHD formalism, where they play the role of the equation of state of the system, completing the evolution equations to allow for the determination of time evolution from the initial conditions.A further advantage is that they have a relatively simple encoding in terms of graphical diagrams, facilitating the construction and programming of the system.These graphs can be determined from a continued fraction expansion of the coupling, which is finite for rational couplings, and its length is called the number of levels.The first of these levels contains massive nodes describing the massive excitations, such as the soliton (corresponding to the kink/anti-kink) and the breathers; subsequent levels involve the so-called magnons that can be considered as auxiliary quasi-particle excitations encoding the charge degrees of freedom of the kink/anti-kink particles.Irrational coupling values can be described by truncating the continued fraction expansion to a level which gives a desired approximation precision.We gave explicit examples of TBA systems for up to 4 magnonic levels.
The graphs we obtained are similar, but not identical, to the well-known sine-Gordon Ysystem [52].Establishing the equivalence of the two systems is an interesting task that we leave for further study.In this work, by computing the central charge, we only checked that our system gives the correct ultraviolet (high-temperature) limit for at most two magnonic levels.
We then applied sine-Gordon thermodynamics to compute the Drude weight for the charge and energy transport.In contrast to the case of the charge Drude weight studied in [54], quantities involving the energy show no fractal structure, which can be understood from the fact that while the reflective scattering of kinks with anti-kinks introduces a random walk component for the transport of the topological charge, this effect is absent for the energy transport.Indeed, a naive expectation would be to find diffusive transport (i.e.vanishing Drude weight) for the charge.In this connection, we note that in the framework of the socalled hybrid semiclassical approach [39], the (partial) inclusion of the non-diagonal form of kink scattering does indeed lead to diffusive effects [40,51,84].The ballistic transport we find strongly suggests the existence of yet unknown conserved quantities that are odd under charge conjugation, paralleling the case of the XXZ spin chain [68].
Our second application was to the full GHD system, where we considered the simple but paradigmatic setting of a bipartitioned initial state consisting of two semi-infinite parts, each in a different thermal equilibrium state.We demonstrated that the usual iterative scheme [63,80] can be used in both the attractive and repulsive regimes of the sine-Gordon model to obtain the energy and charge (and their current) profiles.We gave an example where the discontinuities [81], coming from the effective velocities of the magnons being less than the speed of light, are shown.Furthermore, we also performed an alternative computation of Drude weights using the bipartition protocol with infinitesimally imbalanced initial states.The results of this calculation match the results of the direct TBA calculation to a very high precision, further supporting the self-consistency of the TBA system and the correctness of the methods used.
We also calculated dynamical correlators in thermal states.The shape of the correlators depends strongly on the temperature because of the interplay of the effective velocity and the densities of occupied states.More interestingly, the support of the charge-charge correlators clearly differs from that of the energy-energy correlators.This is related to a novel effect discussed elsewhere [83].
In addition to the results presented here, the sine-Gordon GHD system completed by the TBA system and the associated dressing equations opens the way to study the sine-Gordon model's hydrodynamics at generic coupling values, which offers many potential applications.Beyond the partitioning protocol studied here, more general inhomogeneous setups are expected to lead to further results relevant to condensed matter and cold atom experiments.In particular, the sine-Gordon model is realised in atom chip experiments [23] in which generalised hydrodynamics can be investigated [85,86] due to the ability of phase imprinting and of designing arbitrary inhomogeneous optical potentials [87,88].Using ballistic fluctuation theory [89][90][91], GHD can also access fluctuations and full distribution functions of conserved charges and currents.Additionally, dynamical correlation functions of vertex operators can also be computed.Studying diffusive corrections to the ballistic behaviour and the possibility of superdiffusive transport [92,93] also provide promising avenues of further research, which we plan to address in the near future.At the reflectionless points ξ = 1, 1/2, 1/3, . . ., there are no magnonic degrees of freedom, and the number of breathers is given by n B = ξ −1 − 1.At the same time, the soliton and antisoliton excitations enter separately and symmetrically.In this case, p 0 = α = 1, resulting in a single kernel with the Fourier transform which in real space has the form The decoupled form of the equations was obtained in [49] with the result The above TBA system has the general structure Eq. (26), where the coupling between the degrees of freedom can be encoded by a graph shown in Fig. 11, with each link corresponding to the same kernel (B.2), while the constants appearing in (26) are given in Table 5.

B.1 Decoupling at level 1
For a system with a single magnonic level, the coupling can be written as Table 5: Parameters appearing in the TBA system Eq.( 26) and the dressing equation (47) at reflectionless points.
Excitations Labels Antisoliton S m S cosh(θ )/T + µ/T +1 +1 +1 The decoupling calculations are easiest to work out in Fourier space, where the convolutions become multiplications, and the computation reduces to simple, albeit somewhat tedious, algebraic matrix manipulations [95].These cases include the points where ξ is a positive integer [81] by setting n B = 0.For general couplings, the coupled TBA equations can be written in Fourier space as where the dashed lines delineate the row/column with index n B + 1, Φ η i j = η j Φ i j and L i = log 1 + e −ε i .We also separated the parts depending on the masses and the charges (i.e.temperature and the chemical potential) in the source terms as w i = w m,i + w q,i .We define the following matrix while all other elements in row n B + 1 are zero.It turns out that similarly to M, the matrix M(1 − Φ η ) has a sparse structure, and we define new matrices K (1) and K (2) as M = 1 − K (1) , M(1 − Φ η ) = 1 + K (2) .(B.10) Multiplying (B.5) by M from the left yields ε i − w m,i − w q,i − K (1) i j w q, j = K (1) i j ε j − w m, j + L j + K (2) i j L j .(B.11) Note that the charge part of the driving term w q,i is constant in rapidity space, so i j w q, j = K (1) i j (t = 0)w q, j δ(t = 0).We define the modified source terms w m,i = w m,i , w q,i = w q,i − K (1) i j w q, j , w i = w m,i + w q,i .(B.12) In writing down the general TBA structure, we also exploit that w m is zero for magnons and that w q,S turns out to be always zero.We now consider the cases ν 1 = 2 and ν 1 > 2 as they require separate treatments.

The case ν 1 = 2
This case was already considered in Ref. [50].The relevant matrices are the following: To get to the final form of the system, the last equation for magnon M 2 can be used to rewrite the soliton equation S Putting all together, the resulting K matrix is Figure 13: Graphical representation of the TBA system at one magnon level and ν 1 = 3 and ν 1 > 3. The various links encode kernels as specified in Table 1.
The graph describing this matrix is shown in Fig. 12, while the other ingredients of the TBA system Eq.( 26) are listed in Table 6.In Fig. 12 and all subsequent graphs, we introduced the common convention that filled nodes denote massive excitations, while empty nodes correspond to massless ones (magnons).

The case ν 1 > 2
In this case, the procedure is the same.However, the end result is slightly different in detail: , basis: , (B.17) which is depicted in Fig. 13.The other ingredients appearing in the TBA system at one magnonic level are summarised in Table 6.Table 6: Source terms and coefficients appearing in the TBA system Eq.( 26) and the dressing equation ( 47) for one magnonic level.For

B.2 Decoupling at level 2
For a system with two magnonic levels, the coupling can be written as This case includes repulsive regime couplings ξ = ν 1 + 1 ν 2 by omitting breathers n B = 0. Note that in this case, it is possible that ν 1 = 1, which must be treated separately for ν 2 = 2 and ν 2 > 2. The computation itself is similar to the level 1 case, and the resulting TBA kernels are specified by the graphs in Fig. 14 for ν 1 = 1 and in Fig. 15 for ν 1 ≥ 2. Note that for ν 1 = 1 the soliton gains a self-coupling with an additional negative sign compared to the generic kernel indicated by the loop turned upside down, and additional negative signs appear in the kernels connecting the soliton to the first magnon.
The other ingredients appearing in the TBA system at two magnonic levels are summarised in Table 7.
Table 7: Source terms and coefficients appearing in the TBA system Eq.( 26) and the dressing equation (47) for two magnonic levels.

Figure 1 :
Figure1: The number of species (breathers, one or two solitons, and magnons altogether) appearing in the Bethe Ansatz equations as a function of the coupling strength shows an intricate structure.We adopt the convention of denoting couplings corresponding to zero (reflectionless), one, two, and three magnonic levels with green circles, blue squares, red diamonds and brown stars, respectively.

Figure 7 :
Figure 7: Comparison of the energy and charge Drude weights at different temperatures as a function of the coupling strength.Even though the number of particles has an intricate dependence on the coupling (cf.Fig.1), the TBA system captures the continuous coupling dependence of the energy Drude weight.In contrast, the charge Drude weight shows a fractal/popcorn dependence on the coupling strength.Besides the numerically calculated values, the bottom left panel shows Eq. (54) with a red dashed line, showing excellent agreement at the reflectionless points.In the rightmost figures, we omit individual points for clarity and only display the data using continuous black lines, as the fractal structure is not resolved numerically at these high temperatures.The continuous behaviour of the charge Drude weight (bottom right panel) in the limit T → ∞ is justified by the agreement with the result (55) of the analytical calculations in Sec.4.3, shown as a red dashed line.The analytical high-temperature limit of the energy Drude weight, Eq. (57), plotted as a red dashed line in the top right panel, also shows good agreement with the numerical results.

Figure 8 :
Figure 8: Gradual suppression of the fractal structure of the Drude weight with increasing chemical potential.

2 B 3 Figure 11 :
Figure 11: Graph representation of the decoupled TBA system at reflectionless points.
.6)While the column n B + 1 is zero except for the element M n B +1,n B +1 which is 1, the elements in row n B + 1 depend on ν 1 and n B .In all cases, B +1,n B = − Φ p 1 (t) , (B.8)and for ν 1 = 2 there is another nonzero entry,M n B +1,n B +3 = Φ p 1 (t) , (B.9) Figure12: Graphical representation of the TBA system at one magnon level and ν 1 = 2.The various links encode kernels as specified in Table1.Filled nodes denote massive particles, while empty nodes correspond to massless ones.