# Generalized hydrodynamics of classical integrable field theory: the sinh-Gordon model

[-12mm] KCL-MTH-17-07

Generalized hydrodynamics of classical integrable

[0.2cm] field theory: the sinh-Gordon model

Alvise Bastianello, Benjamin Doyon, Gerard Watts, Takato Yoshimura

SISSA & INFN, via Bonomea 265, 34136 Trieste, Italy

Department of Mathematics, King’s College London, Strand, London WC2R 2LS, U.K.

Using generalized hydrodynamics (GHD), we develop the Euler hydrodynamics of classical integrable field theory. Classical field GHD is based on a known formalism for Gibbs ensembles of classical fields, that resembles the thermodynamic Bethe ansatz of quantum models, which we extend to generalized Gibbs ensembles (GGEs). In general, GHD must take into account both solitonic and radiative modes of classical fields. We observe that the quasi-particle formulation of GHD remains valid for radiative modes, even though these do not display particle-like properties in their precise dynamics. We point out that because of a UV catastrophe similar to that of black body radiation, radiative modes suffer from divergences that restrict the set of finite-average observables; this set is larger for GGEs with higher conserved charges. We concentrate on the sinh-Gordon model, which only has radiative modes, and study transport in the domain-wall initial problem as well as Euler-scale correlations in GGEs. We confirm a variety of exact GHD predictions, including those coming from hydrodynamic projection theory, by comparing with Metropolis numerical evaluations.

December 15, 2017

## 1 Introduction

Integrability has provided the backbone for a wide array of exact results in theoretical physics, in as diverse frameworks as quantum chains and classical fields. Despite its development over many years, however, until recently it remained unknown how to access dynamical quantities out of equilibrium efficiently. What happens to integrable many-body systems in situations where homogeneity and stationarity are broken, and where non-trivial currents exist? This problem has become of high importance in particular in the quantum realm, thanks to the advent of cold-atom experiments [1, 2, 3, 4, 5].

The lack of homogeneity breaks most of the standard structures at the core of many-body integrability, such as the inverse scattering method. Nevertheless, these standard structures may be used in conjunction with the idea of emergence of hydrodynamics. With weak, large-scale inhomogeneity and non-stationarity, one describes a many-body state in terms of fluid cells, that is, mesoscopic regions which are assumed homogeneous, stationary and very large compared to microscopic scales, where entropy is maximized. The corresponding Euler-scale dynamics is a consequence of microscopic conservation laws, and translates to a dynamics for the Lagrange parameters of each fluid cell. Within fluid cells, the usual techniques of integrability can be applied. The theory that implements these ideas in the context of quantum integrability is generalized hydrodynamics (GHD) [6, 7] (see also [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] for recent developments and applications), a hydrodynamic theory based on the observation that in quantum integrable systems, entropy maximization gives rise to generalized Gibbs ensembles (GGEs) [26, 27, 28, 29].

The goal of this paper is to extend these ideas to integrable classical field theory. Focussing on the classical sinh-Gordon model, we confirm that GHD can indeed be applied to classical fields, thus further suggesting its large universality. We explicitly compare GHD predictions with Metropolis simulations of classically fluctuating initial states evolved deterministically with the field’s equations of motion. We study both the partitioning protocol (or domain wall initial condition) [30, 31, 32, 33, 34, 35, 36] (see also [6, 7, 11, 12, 13, 18, 23, 37], where two halves are initialized in different homogeneous GGEs and suddenly joined at one point, and Euler-scale dynamical correlation functions in homogeneous thermal states. Our study provides in particular the first numerical tests for the recent GHD constructions of Euler-scale correlation functions in integrable models [17, 24].

The most powerful formulation of GHD is obtained in the quasi-particle picture. This picture is natural from Bethe ansatz integrability of quantum systems, where quasi-particles relate to Bethe roots: in this formulation, GHD is based on the formalism of the thermodynamic Bethe ansatz (TBA) [38, 39, 40]. The resulting GHD equations have also been observed to apply to certain classical gases such as the hard-rod gas [41, 42, 43, 12] or soliton gases [44, 45, 46, 47]. In these cases, the quasi-particles represent the solitons themselves, or any dynamical objects with equivalent scattering features [16]. The general hydrodynamic principles behind GHD should, however, be applicable somewhat more generally, independently from an a priori underlying quasi-particle concept, such as in integrable field theory.

The theory of Gibbs states of classical integrable fields has already been developed in [48, 49, 50, 51], and is easy to extend to GGEs. Interestingly, it takes a very similar form to that of the quantum TBA. The exact form of GGEs, and therefore of GHD, depends on the precise type of modes admitted by the integrable model. Recall that in quantum models, fermionic and bosonic modes give rise to different free energy functions involved in TBA [39], and related statistical factors appearing, for instance, in Euler-scale correlation functions [17] (although it is sometimes possible to use both formulations to describe the same system [52]). Likewise, it has been observed that classical gases are associated to Boltzmann-type free energy functions [17]. In classical integrable field theory, it is observed in [48, 49, 50, 51] that many models admit two types of modes: solitonic modes and “radiative modes”, both of which are involved in (generalized) thermalization processes. Solitons have clear particle-like behaviour, giving GHD modes with Boltzmann-like free energy function, like classical gases. By contrast, radiative modes do not display any obvious quasi-particle character in their dynamics, yet, as we will verify, at the level of Euler hydrodynamics, they can be accounted for by the same quasi-particle GHD equations, with the appropriate free energy function for radiative modes. In this paper we only study radiative modes, as the classical sinh-Gordon model does not contain soliton modes.

The radiative free energy function is reminiscent of that of the classical prediction for black-body radiation. Because of the “UV catastrophe”, radiative modes make averages of fields containing high enough derivatives diverge. This is characteristic of the roughness of (generalized-)thermally fluctuating classical fields. We will verify that the GHD equations nevertheless describe correctly the averages of observables that are finite. In thermal states of the sinh-Gordon model, the expectation value and correlation functions of every individual conserved density or current is UV divergent. For our studies of the partitioning protocol and of correlation functions in thermal states, we will therefore focus on the trace of the stress-energy tensor, a combination of energy density and pressure that is UV finite, and more generally on vertex operators (suitable exponentials of the field). In GGEs containing the first (spin-3) non-trivial conserved charge, the energy density and current are UV finite, and thus we also study the partitioning protocol between two such GGEs and evaluate exact profiles of energy density and current.

Our study of vertex operators is based on explicit GGE expectation values obtained via the semi-classical limit of an ansatz recently proposed in the quantum case [53, 54, 55]. Therefore, the present work also provides the first numerical test of this ansatz, albeit within the classical framework. Since vertex operators are generically not conserved densities or currents, their correlation functions go slightly beyond the original proposition made in [17]. However, the hydrodynamic projection methods used in [17] can be applied to non-conserved fields (see [24]), and the present work provides the first numerical verification of such concepts.

The paper is organized as follows: Section 2 presents a summary of the hydrodynamic concepts applied to integrable systems, with particular emphasis on the similarities and differences among quantum and classical systems. Section 3 presents numerical results in the classical Sinh Gordon model and comparisons with analytical predictions from GHD. In particular, section 3.2 studies the partitioning protocol, and section 3.3 considers correlation functions at the Euler scale in homogeneous thermal ensembles. Our conclusions are gathered in Section 4. Finally, three appendices contain more technical details. Appendix A considers the classical sinh-Gordon model as the semi-classical limit of its quantum counterpart, while appendix B mainly supports Section 3.3 providing the correlation functions on the Eulerian scale in the non-interacting limit. Finally, appendix C contains the numerical methods needed to solve the GHD and to directly simulate the model.

## 2 GHD of classical fields

GHD, and hydrodynamics in general, is based on the assumption of local entropy maximization. It is therefore naturally associated with statistical distributions, such as Gibbs and generalized Gibbs ensembles. Paralleling works done in the quantum case, the natural context in which GHD could be applied to classical fields is not that of the evolution of a single field configuration, but instead that of the deterministic evolution of statistical distributions thereof. This has been considered for instance for anharmonic chains recently [56], where Gibbs ensembles give probability distributions for initial chain configurations, and conventional hydrodynamics successfully describes the evolution of local averages under deterministic dynamics.

Although most of the studies of classical integrable field theory concentrate on exact solutions of the field equations from single field configuration, such as solitonic solutions, there has been important work on statistical distributions of integrable fields as well. There are two categories of results: *i)* those dealing with distributions of solitons [44, 45, 46, 47], and *ii)* those dealing with Gibbs ensembles of classical fields [48, 49, 50, 51]. Soliton gases, at the Euler scale, also satisfy hydrodynamic equations, which were in fact observed to be those of GHD in [16]. Hence one may expect the GHD of soliton gases to play a role in the hydrodynamics of classical field theory with solitonic modes. However, Gibbs ensembles of classical fields display, in general, both solitonic and radiative modes [48, 49, 50, 51].
Below we extend these to generalized Gibbs ensembles (GGEs) (this is a very simple extension). From this we then obtain, using the same hydrodynamic principles as those recalled in [6, 7], GHD for classical fields. The solitonic part agrees with the hydrodynamics of soliton gases, but the radiative part gives important extra contributions.

### 2.1 GGEs

We consider a field theory (relativistic or Galilean) with Hamiltonian , depending on canonically conjugate fields and . For the purpose of the general discussion, these may take values in , in or in any other target space.

This theory is assumed to be integrable, with a space of conserved charges in which lies, spanned by some basis for . In the conventional view on integrability, this would be taken as a basis of local conserved charges: charges that can be written as where is a local functional of and at (that is, involving these fields or their derivatives and products thereof at the point ). The modern viewpoint also admits the so-called quasi-local conserved charges. These have densities with “finite support” (usually with exponentially decaying tails away from ). They have been fully understood in the quantum XXZ spin chain [57, 58, 59], and, completing with respect to an appropriate inner product on the space of conserved charges, they generate the Hilbert space of pseudo-local charges, rigorously shown in quantum lattices to be involved in generalized thermalization [60]. The general picture is similar in the case of quantum field theory, although the understanding of non-local conserved charges has not reached the same level of maturity [61, 62, 63, 64, 65, 66, 67, 68]. It is natural to expect that quasi-local charges should also be involved in GGEs arising from non-equilibrium classical field theory. For the present purpose, we simply follow the standard arguments of GHD [6, 7] and assume appropriate completeness of in order to obtain GHD equations.

We are interested in generalized Gibbs ensembles. These are statistical distributions of field configurations with averages described formally as^{1}

(1) |

for any (quasi-)local functional of and at . Gibbs ensembles at inverse temperature , with instead of in (1), were studied in the sine-Gordon model [48, 49, 50] and sinh-Gordon model [51] using the classical inverse scattering method, and in some generality using the semi-classical limit of the quantum Thermodynamic Bethe ansatz [69]. These studies gave rise to expressions for the free energy that closely resemble those found in the thermodynamic Bethe ansatz of quantum integrable models [38, 39]. Extending to generalized Gibbs ensembles, this in turns allows for the evaluation of averages of (quasi-)local conserved densities in classical GGEs (here and below, a field written without explicit space-time position is assumed to be at the origin).

Under the evolution determined by the equations of motion of the field theory in the usual fashion, conserved densities satisfy the conservation equation

(2) |

The currents are also (quasi-)local functionals of and at . Their averages in GGEs play a fundamental role in the hydrodynamic description. These do not immediately follow from the free energy. However, the point of view in which the classical field theory is the semi-classical limit of a quantum field theory – which is certainly valid in the sinh-Gordon and sine-Gordon models for instance – allows us to take the fundamental results of [6, 7] in order to arrive directly at expressions for the averages of conserved currents in classical field theory.

We now state the main results. For the purpose of clarifying the structure, we first put the results within a general framework, abstracting the structure of the formulae found in [48, 49, 50, 51] for the sine-Gordon and sinh-Gordon models. As per the above comments, these are natural extensions to GGEs, and to averages of currents, of the formulae found there.

In the present viewpoint, the inverse-scattering solutions of the field theory involves a set of possible quasi-particles. Here we use the “quasi-particle” terminology in a loose way: each quasi-particle may be of solitonic or radiative type, the former representing isolated poles in the scattering data, the latter singularities that are part of a continuum (such as a branch cut). Quasi-particles are also characterized by rapidities lying in , which in the Galilean case can be thought of instead as velocities. The doublet can be seen as a “spectral parameter” for the quasi-particle, and the spectral space. An inverse scattering solution is fully determined by giving a set of spectral parameters of solitonic type, and densities over spectral parameters of radiative type. This has the physical meaning that the solution contains these quasi-particles (solitons and radiative modes) with these rapidities and densities. Each conserved charge is fully characterized by a spectral function . This is in the sense that, when evaluated in a field configuration corresponding to the set and the densities , it takes the value

(3) |

where the sum over is only over quasi-particle type of radiative type, and the spectral parameters are all of solitonic type. For instance, a model composed of independent free fields contains radiative modes and no solitonic modes, and the densities (for ) are simply related to the Fourier transforms of the fields. By contrast, the sine-Gordon model contains two solitonic modes (the soliton and anti-soliton) and one radiative mode. See [69, 70, 71, 72, 73] for details of the inverse scattering method.

Recall that a GGE is characterized by a set of generalized inverse temperatures . It was found in [48, 49, 50, 51] that, as in the quantum case, the main object for the thermodynamics of classical field theory is the pseudo-energy , which depends on a spectral parameter . This solves the following integral equation [48, 49, 50, 51]:

(4) |

where the source term is determined by the generalized inverse temperatures as

(5) |

and where is the free energy function, which is a characteristic of the type of mode that the quasi-particle represents:

(6) |

For the sake of comparison, recall that the pseudo-energy of the quantum TBA is of a similar form, where the free energy function is instead

(7) |

In (4), the quantity , which we will also denote by for and , is the “differential scattering phase”, which characterizes the interactions between the quasi-particles (we assume it to be symmetric). See [48, 49, 50, 51] for its explicit form in the sine-Gordon and sinh-Gordon models (the latter is recalled in the next section).

As usual, the model is further specified by two spectral functions: the momentum and the energy (e.g. in relativistic models and where is the mass of the quasi-particle ). Averages of local conserved densities are evaluated by differentiating the free energy

(8) |

with respect to (here and below the prime means rapidity derivative, e.g. ). Averages of currents are similarly obtained by differentiating [6]

(9) |

The results are

(10) |

and

(11) |

where we use the notation

(12) |

In these expressions, we have introduced the usual TBA quantities and operations: the occupation function

(13) |

the dressing operation

(14) |

the quasi-particle density

(15) |

and the effective velocity

(16) |

Again, (13) is to be contrasted with the expressions in the quantum case,

(17) |

The second expression on the right-hand side in (10) is to be compared with the expression (3) of the conserved charge on a single field configuration: we see that the solitonic modes are now also described using densities, representing the fact that we have a thermodynamically large number of solitons. We do not know of an equivalent expression to in (11) for the conserved current on a single field configuration in terms of inverse scattering data.

### 2.2 Ghd

Given that the structure of GGEs in classical field theory is the same as that for quantum models, up to a modification of the free energy function, it follows that all GHD results can immediately be adapted to the hydrodynamic of classical fields, under the same hydrodynamic assumption. That is, consider local observables in states with weak, large-scale inhomogeneities:

(18) |

Here the time-evolved field is obtained from the equations of motion of the field theory in the usual fashion. If vary non-trivially only on very large scales, the main assumption is that

(19) |

with space-time dependent GGEs,

(20) |

The large-scale conservation laws give the GHD equations. In terms of the fluid variable this boils down to [6, 7]

(21) |

Combined with the description of GGEs in the previous subsection, this is the Euler hydrodynamic equation for classical field theory.

We can now immediately apply the various results from GHD. Below, we will make use of two sets of results: those pertaining to the partitioning protocol [6, 7], and those relating to Euler-scale correlation functions in GGEs [17, 19, 24]. The former give averages in non-equilibrium states where non-zero currents exist, while for the latter we will restrict ourselves to dynamical correlations in homogeneous, stationary states (see [24] for the case of inhomogeneous, non-stationary states). We briefly review the main results.

In the partitioning protocol, the initial state has probability measure

(22) |

This is not of weak variation, as there is a sharp jump at the origin. However, after a small relaxation time, it has been observed both in quantum spin chains [7] and in classical gases [12] that the generalized hydrodynamic description applies, with the corresponding hydrodynamic initial state. We thus set (where is Heavyside’s step function), which is composed of a GGE on the left-hand side, , and a GGE on the right-hand side, , and we solve (21) with this initial condition. The solution depends only on the ray . It is obtained by solving the following “self-consistent” set of integral equations, for the functions and :

(23) |

Here, is determined by (16) where the dressing is with respect to the occupation function . This set of equations was found to be solvable rather efficiently by iteration [6, 7].

GHD also gave rise [17, 19, 24] to formulae for connected correlation functions

in homogeneous, stationary GGEs at the Euler scale. As we numerically observe below, Euler-scale correlation functions are more precisely obtained by averaging over fluid cells. We expect there to be many ways of doing such averaging. Consider for instance a neighbourhood around the point , which can be thought of as a disk of radius that grows with under the condition . Denote its area by . Then, GHD predicts

(24) |

Here is the statistical factor, which equals

(25) |

This is again to be compared with the quantum case,

(26) |

The functions are obtained from the hydrodynamic projection

(27) |

which is valid for any local observable. In the case of charge densities and currents, the explicit expressions are and .

Of course the integral in (24) can be performed and is supported, thanks to the delta function, on the values of for which . Integrating over space , the factor disappears and we obtain

(28) |

where we keep the time-domain fluid-cell average, with an interval of radius centred on . See [17, 24] for more details.

Remark. As mentioned above, the GHD of classical fields can be seen as composed, in general, of two types of “quasi-particles”: solitonic modes and radiative modes. All aspects of GHD that pertain to the solitonic modes – such as the solitonic part of the formula for effective velocity and the statistical factors in correlation functions – agree exactly with the corresponding aspects of the GHD for gases of solitons. That is, there is a part of the hydrodynamics which arises from the gas of the field theory’s own solitons. However, this is not enough to describe correctly the hydrodynamics, as the radiative modes are also necessary. We also emphasize that the statistical factor of solitonic modes agrees with the statistical factor for particle gases, such as the hard-rod gas [17]. Thus solitons truly behave as particles, while radiative modes do not.

## 3 The sinh-Gordon model

The results presented above follow from *i)* generalizing the results of [51, 69], for the sine-Gordon and sinh-Gordon models, to arbitrary models in arbitrary GGEs, and *ii)* combining with the results of [6, 7] in order to get the hydrodynamics. The generalization to GGEs is a straightforward exercise, and the abstraction to arbitrary models is a natural guess. The expressions for the current averages in (11), and the effective velocity (16), were never written down before in classical field theory. However, they are again natural from the fact that GGEs, in the quasi-particle formulation, have a very similar structure in classical and quantum field theories. For completeness we present in Appendix A a full derivation of the hydrodynamics via the semi-classical limit of the quantum sinh-Gordon model, following the ideas of [69]. Since in the quantum sinh-Gordon model the corresponding equations, and in particular the effective velocity (16), were derived in [6], this gives a full derivation of all necessary ingredients for the hydrodynamics of the classical sinh-Gordon model (it would be interesting to have an independent derivation of the effective velocity fully within classical field theory, which we leave for later works).

An advantage of dealing with a classical theory, as compared with a quantum model, is the relative simplicity of testing it through numerical simulations. This is especially true for continuous models, where the DMRG methods [74] are not efficient. In this section, after specializing the hydrodynamics to the classical sinh-Gordon model, we test its predictions against direct numerical simulations, whose details are left to Appendix C.

### 3.1 GGE of the sinh-Gordon model and UV finiteness

The sinh-Gordon model (ShG) is described in terms of a single scalar field whose dynamics is ruled by the following Lagrangian:

(29) |

The model is shown to be classically integrable by means of the inverse scattering method [70, 71, 72, 73, 69] and it is a well known fact that its integrability survives under quantization [75, 76]. At the classical level the spectrum of ShG possesses only one purely radiative species [69], making its TBA and the subsequent hydrodynamics rather simple. The resulting pseudo-energy equation [51, 69] can be brought to the form (4), in the case of radiative mode as per (6), with , or equivalently

(30) |

with

(31) |

The integral in (30) is a Cauchy principal value integral for the singularity at . In [69], (30), (31) are derived both through the inverse scattering method and the semi-classical limit of the quantum TBA^{2}

It is well known that thermal states in classical electromagnetism suffer from UV catastrophes, intimately connected to Planck’s original idea of quantizing oscillation modes in order to explain the black body radiation spectrum. The same problem appears in thermal ensembles and GGEs of other field theories such as the sinh-Gordon model. The problem is that in such states fields can be very “rough”. Essentially, a path integral of the form (1), if it involves local conserved charges only up to some finite spin, does not regularize the field enough to guarantee the existence of averages of its large-order derivatives. In the TBA formulae above, this is clearly seen as a divergence of averages of conserved densities of high enough spin, which occur at large rapidities due to radiative modes. Indeed, consider (10) and (11) in the sinh-Gordon model. At large rapidities, is bounded. However, where is the source term in (4), while . Thus, convergence of the average of the density or current corresponding to is guaranteed only if is integrable. In a thermal ensemble, (where is the inverse temperature), and therefore no local conserved density or current has finite average (except for those that are zero by parity symmetry) – recall that in the sinh-Gordon model, local conserved densities are those with with odd. Nevertheless, it is a simple matter to see from the TBA formulae that the right combination of energy density and momentum current is finite: this is the trace of the stress-energy tensor , which does not contain field derivatives. Clearly, averages of higher-spin densities and currents become finite as we increase the spin of the conserved charges involved in the GGE.

Vertex operators also have finite thermal averages, as they do not contain field derivatives.
Since GHD provides the full GGE state at every space-time point – not just the averages of conserved densities and currents – we can combine GHD information together with homogeneous GGE results to study these quantities. For instance, following [69], we could access their expectation values trough the LeClair-Mussardo expansion [77], but this provides a small excitation-density description, which is not best suited for the present purpose. However, an ansatz for the ratios of expectation values of vertex operators in the *quantum* model, valid for arbitrary excitation density, exists. It has been firstly proposed for thermal ensembles [53, 54] and recently extended to arbitrary GGEs [55]. Taking the classical limit, we obtain formulae for expectation values in the classical sinh-Gordon model; this is a new result within the classical realm. The derivation is reported in Appendix A.

### 3.2 The partitioning protocol

*Analytical expectation values of the trace of the stress energy tensor and of vertex operators in a homogeneous thermal ensemble at inverse temperature are compared with the free model. With a suitable rescaling of the fields, the expectation values depend only on the combination , thus we choose and vary . First panel on the left: the traces of the stress energy tensor for the interacting and free cases are compared. Second and third panel: the same is done for vertex operators. Note that passing from the interacting to the free theory, the trace of the stress energy tensor changes as function of the field. The difference between interacting and free is enhanced when comparing expectation values of fixed functions of the fields.*

In the partitioning protocol, the initial condition is set as per (22). For definiteness, we set to be the energy density, the energy current, equal to the momentum density, and the density of the first non-trivial (spin-3) conserved charge (see (107) for its explicit expression, and (110) for its expression in the light-cone lattice regularization of the model). We set for all . We use (4) and (13) in order to determine and , where and are the one-particle eigenvalues of the energy and of the spin-3 charge, respectively.

We study two types of observables. The energy density and current , have averages predicted by (10-11) along with (23). In particular, the trace of the stress-energy tensor – the difference between the energy density and the momentum current – is given by

(32) |

Note that given a relativistically invariant Lagrangian written in terms of a single real scalar field , the canonical stress energy tensor is defined as

(33) |

with the Minkowski two dimensional metric. Thus is (twice) the ShG potential

(34) |

We also study generalizations of this, the vertex operators . Based on Ref. [53, 54, 55], a recursive relation for the GGE averages of vertex operators can be obtained (for details see Appendix A)

(35) |

where

(36) |

Taking to be an integer (and using the fact that if ), we obtain explicit formulae for .

We consider two cases of the partitioning protocol: the purely thermal case, with non-zero and , and the GGE case, the situation with higher conserved charges in the initial left and right states, and all non-zero. As explained above, in the purely thermal case, because of a UV catastrophe similar to that of black body radiation, the expectation values of energy and momentum densities (as well as all the individual local charge densities and currents) are divergent, but the averages of the trace of the stress-energy tensor – the difference between the energy density and the momentum current – and of the vertex operators are finite, as these do not contain field derivatives. Finiteness of the GGE average of is also clear from the integrability of the integrand in (32). Thus in the thermal case we concentrate on these observables. In the GGE case, we study the energy density and current independently.

*Profiles of the trace of the stress energy tensor and of the first vertex operators as functions of the ray . Parameters: , , inverse temperatures of the initial ensembles and , the number of samples of the Metropolis is 84000. An additional average along the ray from up to smears out the statistical oscillations. In terms of the dimensionless coupling , we are in the strongly interacting regime (see Figure 1).*

We compare the analytic predictions with direct numerical simulations of the out-of-equilibrium protocol, see Appendix C for details of the numerical methods. In order to illustrate the effect of the interaction and determine the values of which significantly affect the averages, in Figure 1 we plot the analytical thermal averages of the stress energy tensor and of the vertex operators both in the interacting and the free () purely thermal Gibbs ensemble at the same temperature.

The results of the comparison between GHD analytical predictions and numerical simulations in the thermal case are presented in Figures 2. We find excellent agreement, with the measure of difference between the numerical and analytical values defined by the ratio (where is the numerical data and is the analytical value) given by , and in the cases of Figures 2a, 2b and 2c respectively. Our data is precise enough to show a clear departure from the free theory result, see Figure 3. In Figure 4 we do the comparison for the energy density and energy current in the GGE case. Again very good agreement is found with the measure of error being and for the charge and the current respectively.

It must be stressed that, despite various checks in known limiting cases (see [55]), the formula for the vertex operators (35-36) had never been numerically verified before. Figure 2 constitutes both a verification of the validity of generalized hydrodynamics for the sinh-Gordon model *and* a numerical verification of the classical limit of the vertex operator ansatz.

*(a) and (b) are numerically computed (red line) and compared with the GHD predictions (black line). The numerical results are obtained by averaging over the range . The parameters are with the coupling to the higher-spin charge varying, , and the lattice spacing in the light-cone discretization is*

### 3.3 Euler scale correlation functions

Equation (24) gives a prediction for the behaviour of the connected two-point correlators of local observables at the Euler scale: at long times and large distances, after fluid-cell averaging. Such expressions in integrable systems have never been directly compared with numerical simulations and here we present the first study in this direction. We consider correlations in homogeneous thermal ensemble. The UV catastrophe still prevents us from computing correlation functions of arbitrary charges and currents, however one can check from the TBA formulae that a certain number of stress-energy tensor correlation functions do stay finite, such as and the difference , and similarly certain correlation functions involving vertex operator are finite. In the following we focus on two correlation functions: the trace two-point function , and the symmetrized vertex operator two-point function , because these offer the best chance of numerical precision.

We first remark that, without fluid-cell averaging, one may expect persisting oscillations, of the order of the hydrodynamic value itself, at large space-time variables . Such oscillations are beyond the validity of the hydrodynamics, which becomes predictive only when fluid cell averages are considered. This is clear at the free-field point : the large-time limit with and fixed, of the correlation function rescaled by displays non-vanishing oscillations around the predicted hydrodynamic value, see Appendix B. In the interacting case, the presence of such persisting oscillations is much harder to assess. In Figure 5 we study the scaled correlation functions

(37) |

and

(38) |

on rays at various times. Relatively large oscillations do arise at all times we were able to reach, although it is inconclusive whether these actually persist asymptotically.

In order to integrate out these oscillations most efficiently, we study the time-averaged and space-integrated scaled correlation functions (24). More precisely, on the one hand time averages along rays directly give the Euler-scale correlation function,

(39) |

We have numerically observed, and analytically calculated in the free case, that it is not necessary, after such time averaging, to perform any mesoscopic space averaging in order to recover the Euler-scale hydrodynamic result: the large time limit as in (39) indeed exists and gives the right-hand side.

On the other hand, space integrations are obtained as (28). In this case, we have analytically observed at the free-field point that, for certain observables, mesoscopic time averaging (as in (28)) is necessary in order to recover the predicted hydrodynamic result (see Appendix B). That is, at least for correlators in the generic form , undamped time-oscillating terms survive the space integration, but are cancelled by a subsequent fluid-cell time average. In the interacting case, again our numerics does not allow us to reach unambiguous conclusions, although we observe long lived oscillations of the space integrated correlators around the GHD prediction within the numerically accessible time scale. It turns out, however, that mesoscopic time averaging is not necessary for the correlation functions and that we have chosen; this is seen analytically in the free case and numerically observed in the interacting case, and simplifies our numerical checks.

*The predictions (black curves) for the time average of two-point functions of (a) the trace of the stress-energy tensor, (40), and (b) vertex operators, (39) with (42), are compared to the numerical simulation (red curves). In the case (a) the free result of Appendix B (dotted curve) is plotted for comparison. In order to improve the convergence, short times are avoided: the time average is taken in the interval , rather than as in eq. (40), with and . The same parameters as those of Figure 5 are used.*

For the trace of the stress-energy tensor, time averaging along rays gives

(40) |

while space integration gives

(41) |

We refer to Appendix B for a direct analytical check of (40) and (41) in the free case, together with an analysis of the leading corrections. For the vertex operators, the hydrodynamic projection functions in (39) are calculated using the GGE one-point averages along with (27). We show in Appendix A that the function is given by

(42) |

where , and where is the solution of the following linear integral equation:

(43) |

Then one directly applies (39) and (28). Again we refer to Appendix B for a direct analytical check of vertex operator correlations functions in the free case.

In Figures 6 and 7 the time-averaged and space-integrated correlators, respectively, are considered. We see very good agreement between the numerical simulation and the hydrodynamic prediction. In Figure 8 we study space-integrated correlation functions of vertex operator without symmetrizing. In this case, as discussed above, we observe long-lived oscillations around the hydrodynamic value.

*The predictions (black curves) for the space integral of two-point functions of (a) the trace of the stress-energy tensor, (41), and (b) vertex operators, (28) with (42), are compared to the numerical simulation (red curves) as functions of time. At large times, the oscillating ripples are due to the Metropolis noise and are observed to reduce on increasing the number of samples. The same parameters as those of Figure 5 are used.*

## 4 Conclusions

In this paper we studied the hydrodynamics of the classical sinh-Gordon model. In the set-ups we studied, the initial state is distributed according to (generalized) Gibbs ensembles, and this distribution is deterministically evolved following the classical field equation. We studied both one-point averages in the non-equilibrium partitioning protocol, and dynamical two-point functions in homogeneous states. We compared predictions from the new theory of GHD with results from a direct numerical simulation of the protocol. We found excellent agreement in all the cases we studied.

This gives direct evidence that GHD, first developed in the context of quantum chains and quantum field theory and later observed to be applicable to classical integrable soliton-like gases as well, is also applicable to classical field theory, thus further confirming its wide applicability. Because of the UV catastrophe of classical field theory, not all local observables have finite averages in thermal states or GGEs. Nevertheless, we verified that for those whose averages are finite, GHD predictions are correct.

The TBA formulation of GGEs in classical field theory generically involves two types of modes: solitonic modes and radiative modes. The latter are responsible for the UV catastrophe of classical fields. Similarly, the GHD of classical integrable field theory is that for a gas of solitons and radiation. Solitonic modes have a natural classical-particle interpretation, and the part of GHD related to such modes is indeed identical to the GHD of classical soliton gases. Radiative modes, on the other hand, do not have a clear particle interpretation, and were never studied within GHD before. Nevertheless the quasi-particle formulation of GHD still holds, with the appropriate radiative free energy function. In the classical sinh-Gordon model, only radiative modes occur, hence our results give a verification of GHD in this new sector of the theory.

Our results also constitute the first direct verification of the recent GHD correlation function formula. Crucially, we analysed not only conserved fields – specifically the trace of the stress-energy tensor, which is UV finite – but also vertex operators, local fields which are not part of conservation laws. GHD two-point functions for such fields involve the hydrodynamic projection theory, and our results confirm its validity. In addition, our results also provide the first numerical verification of GGE one-point functions formulae for vertex operators recently proposed.

In our study of correlation functions, we observed that GHD predictions are valid up to oscillations which may not subside at large times. We have observed this explicitly in the free-field case, where in certain cases, oscillations both in space and time, of the same order as that of the hydrodynamic predictions, persist indefinitely. For certain observables, appropriate fluid cell averaging appears to be essential in order to recover the Euler-scale correlation predictions from hydrodynamics.

It would be interesting to analyze the approach, at large times, to the GHD predictions. We have explicitly evaluated (and numerically verified) large-time power-law correction terms in the free-field case, but preliminary numerics suggest that in the interacting case, some correction terms are instead exponential. It would also be interesting to extend this analysis to the non-relativistic limit of the sinh-Gordon model, namely the well known nonlinear Schrödinger equation, and to other models, such as the sine-Gordon model, which possesses solitonic modes, and the classical Toda chain.

## Acknowledgments

AB is indebted to Andrea De Luca for useful discussions and to King’s College London for hospitality in the framework of the Erasmus Plus traineeship mobility programme. BD is grateful to H. Spohn for discussions, and thanks the Institut d’Études Scientifiques de Cargèse, France and the Perimeter Institute, Waterloo, Canada for hospitality during completion of this work. BD also thanks the Centre for Non-Equilibrium Science (CNES). GW would like to thank N. Gromov for discussions on conserved quantities. TY acknowledges support from the Takenaka Scholarship Foundation.

## Appendix A Semi-classical limit of the sinh-Gordon model

Even though the (generalized) TBA for classical models can be directly accessed by mean of the inverse scattering approach [78], the rather technical derivation appears to be more cumbersome compared with the quantum counterpart, which probably is more familiar to the reader. This appendix provides a self-consistent derivation of the TBA and hydrodynamics for the classical sinh-Gordon, regarding the latter as the classical limit of its quantum version. Concerning the TBA, such a route was already considered in [69], together with the classical limit of the LeClair-Mussardo formula [77] leading to a small excitation-density expansion for one-point functions. Here we resume such a program extending the derivation to the whole hydrodynamics, then we present a close expression for the expectation values of vertex operators obtained through the semi-classical limit of the ansatz presented in [69]. While in agreement with the LeClair-Mussardo formula, this last result is not constrained to small excitation densities and is feasible to be applied to arbitrary GGEs.

The quantum sinh-Gordon is arguably one of the simplest, yet interacting, quantum field theories. Its spectrum consists of only one type of particle [75], leading to a simpler version of the TBA equation (4)

(44) |

where we explicitly use the subscript “q” as a remind we are dealing now with the *quantum* model; in the absence of any subscript we will refer to classical quantities. Differently from [69], along the semi-classical limit procedure we will use the fermionic formulation of the quantum ShG rather than the bosonic one: even though the two descriptions are equivalent [69], the quantum hydrodynamics is formulated in the fermionic basis making the latter a preferred choice. Referring to the parameter choice of the (quantum) Lagrangian (29), the quantum kernel is

(45) |

The physical mass of the particles is related to the bare mass through the relation [79]

(46) |

As already pointed out in [69], the semi-classical limit of the ShG model can be obtained through a combination of a small coupling and high temperature limit. More precisely we consider the rescaling

(47) |

with a generic observable functional of the field and the Lagrange multipliers of the GGE. Then, in the limit, the quantum expectation value reduces to the expectation value of the classical observable on a GGE having as Lagrange multipliers. Such a claim is easily checked in the simplest case of a thermal state: using the path integral representation of the latter we have

(48) |

where the euclidean time runs on a ring of length . In the high temperature limit the integral over the compact dimension can be approximated with the integrand at times the length of the ring. In the same limit the measure collapses to where the fields are now restricted to .

(49) |

A final change of variable in the path integral makes explicit the appearance of the classical expectation value of the observable on a classical thermal ensemble of inverse temperature . A similar reasoning can be extended to more general GGEs as well to the dynamics, confirming (47) as the correct scaling to obtain the desired semi-classical limit. Notice that, in such a limit, the physical quantum mass simply reduces to the bare mass.
Concerning the TBA equation (44), the rescaling of the Lagrange multipliers is simply translated in . A finite integral equation in the limit is obtained provided we redefine , with the classical effective energy (this rescaling differs from that of [69], where the *bosonic* rather than the *fermionic* quantum effective energy is used).
Through straightforward manipulations, at the leading order in the following integral equation is obtained

(50) |

The last step involves an integration by parts that leads to the final independent integral equation

(51) |

where stands for the principal value prescription, natural emerging passing from (50) to (51): this is the sought-after TBA equation for the classical ShG. We should warn the reader about a difference in the sign in front of the integral between the above expression and the same reported in Ref. [69], where a typo must have occurred. This concludes the pure TBA calculation and we can now pass to the hydrodynamics: rather than taking the differential scattering phase for the classical ShG from (51) and plug it into the hydrodynamic formulas of Section 2.2, we directly take the limit of the GHD for quantum integrable systems as a further check of the correctness of its classical version. Through calculations close to those that provided the classical TBA, it is possible to show that the quantum dressing operation (14) (here specialized to the ShG case)

(52) |

Reduces to the sought definition of the classical dressing

(53) |

provided we set while taking the limit. In particular, this implies that the effective velocity (16) of the quantum model converges, in the semi-classical limit, to the definition of the effective velocity in the classical case. Finally, from the definition of the occupation function in the classical (13) and quantum (17) cases and the rescaling of the effective energies we have, at the leading order in , . Replacing this last piece of information in the quantum version of the GHD (21) we readily obtain its classical counterpart. Alternatively, we could have noticed that in the semi-classical limit the expectation values of the quantum conserved charges and their currents go, respectively, to the classical expectation values of charges and currents as written in eq. (10-11) (apart from an overall and inessential factor) and impose the conservation laws.

### a.1 The expectation value of the vertex operators

This part is devoted to the derivation of the equations (35-36) reported in the main text, through the semi-classical limit of the ansatz presented in Ref. [55]. Following the latter, the expectation values of vertex operators in a GGE for the quantum model satisfy

(54) |

with

(55) |

(56) |

As it is evident from the rescaling (47), the quantum expectation value of the vertex operators goes, in the semi-classical limit, directly to the classical expectation value of the latter (the rescaling of the coupling and of the field in the observable cancel each others). Thus the classical limit of the r.h.s. of (54) will directly produce the ratio of the expectation values of the vertex operators in the classical theory. The first step to take the limit is to consider up to first order in . Using that after the rescaling is first order in , we eventually need (in the distribution sense)

(57) | |||

### a.2 The Hydrodynamic projection of the vertex operators

Having access to the exact value of the expectation value of vertex operators on arbitrary GGEs makes the exact computation of the hydrodynamic projection (27) feasible, paving the way to correlation functions on an Euler scale. Specializing (27) to the ShG model and to the vertex operators, we are interested in computing

(58) |

where for simplicity and is the Lagrange multiplier associated with the charge whose eigenvalue is . Differentiating both sides of eq. (35) with respect to , we readily obtain the recurrence relation

(59) |

That has the obvious solution

(60) |

can now be extracted from the integral equation (36). To carry on the calculation, it is convenient to introduce an operatorial representation of the integration kernel. In particular, we define the linear operators and as follows:

(61) |

for any given test function . In this language, the linear equation satisfied by (36) is written as

(62) |

where we introduce . The formal solution is immediately written

(63) |

and the derivative with respect to is easily performed, using the fact that and the TBA equations (30)

(64) |

As last technical step, we consider the integral appearing in eq. (60)