Exact out-of-equilibrium steady states in the semiclassical limit of the interacting Bose gas

We study the out-of-equilibrium properties of a classical integrable non-relativistic theory, with a time evolution initially prepared with a finite energy density in the thermodynamic limit. The theory considered here is the Non-Linear Schrodinger equation which describes the dynamics of the one-dimensional interacting Bose gas in the regime of high occupation numbers. The main emphasis is on the determination of the late-time Generalised Gibbs Ensemble (GGE), which can be efficiently semi-numerically computed on arbitrary initial states, completely solving the famous quench problem in the classical regime. We take advantage of known results in the quantum model and the semiclassical limit to achieve new exact results for the momenta of the density operator on arbitrary GGEs, which we successfully compare with ab-initio numerical simulations. Furthermore, we determine the whole probability distribution of the density operator (full counting statistics), whose exact expression is still out of reach in the quantum model.


Introduction
The out-of-equilibrium physics is a fascinating subject which is currently at the center of an intense research activity: this field has witnessed a striking progress due to both new experimental techniques, mostly available in the world of cold atoms, and an impressive development of analytic and numerical methods (see, for instance, [1]). These advances have opened the door to explore novel phenomena and to answer fundamental inquiries in quantum mechanics and statistical physics, in particular about relaxation and equilibration in many-body systems. One of the most crucial achievements of the last years has been the possibility of experimentally realising almost-perfectly isolated quantum systems and, at the same time, tuning their interactions and the effective dimension so that one can investigate the quantum quench [2] in the lab, as done indeed in several experimental setups [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Arguably, the quantum quench is among the simplest out-of-equilibrium protocols: after the system is prepared in a given initial state, which could be pure or mixed, it is let evolve with a non-trivial time-independent HamiltonianĤ. As discussed below, one of the key topics of the out-of-equilibrium physics concerns with the nature of the late-time steady state which follows a quench, in particular the possibility to predict the asymptotic stationary values of the various observables of the system. In this paper we are going to address the following question: can we characterise the steady state which emerges in the infinite time limit of a classical integrable field theory? As we will show, the answer to this question proves to be remarkably rich and instructive, with several far-reaching consequences for what concerns our theoretical understanding of the outof-equilibrium phenomena both in quantum and classical realms. It is worth explaining why.
The main idea. We take advantage of the interplay between classical and quantum mechanics to further deepen our understanding in both realms. Indeed, as we extensively discuss below, the quantum world has been thoroughly investigated and several important results have been achieved: surprisingly, in several instances tackling directly the quantum system revealed to be simpler than the classical case. To this hand, quantum results can be carried to the classical realm through proper semiclassical limits, achieving new insights. On the other hand, many unsolved questions in the quantum case beg to be addressed: one among all is the determination of the steady-state after a quantum quench, which still lacks a general method for generic initial states (see Ref. [19] and references therein). As we will show, this problem can be instead efficiently semi-analytically solved in the classical case, providing a complete solution of the quantum quench problem, within the semiclassical approximation.
Classical non-relativistic integrable models can be accessed thanks to the attractive features of quantum relativistic integrable field theories (see, for instance [20] and references therein), in particular the possibility to obtain for these theories exact information on their elastic and factorized S-matrix amplitudes [21][22][23], the spectrum of their excitations, the thermodynamics [24,25], the exact expressions of the matrix elements of local operators [26][27][28], together with their correlation functions at zero temperature [29,30] and at finite temperature [31][32][33][34].
Since in addition to genuine coupling constants, any quantum field theory has inherently two important parameters -the speed of light c and the Planck constant ħ h -it should be possible to use the solvability of the integrable quantum field theories to get exact results for all those models which emerge by playing with these two parameters. This will be extensively used hereafter.
Three integrable models. It is useful to introduce the three integrable models which will accompany us in the rest of the paper.
• The first is the relativistic integrable quantum field theory, called the sinh-Gordon (ShG) model [22,25,28], with Hamiltonian given by , where Φ = Φ(t, x) is a real quantum scalar field which satisfies the canonical equal-time commutation relations Of course, replacing commutators with Poisson brackets, [, ] → iħ h{, }, and regarding Φ(t, x) as a classical field, one has the classical ShG model.
• The second is the integrable non-relativistic quantum field theory, called the Lieb-Liniger (LL) model [35][36][37], with Hamiltonian given bŷ where the complex Bose field ψ(t, x) satisfies the canonical commutation relations In the Hamiltonian (2), the interaction is rescaled in order to attain a meaningful semiclassical limit, as we will extensively discuss: while sending ħ h → 0, κ must be kept constant. In the following, we consider only the case κ > 0, i.e. the repulsive regime of the Lieb-Liniger model. The Fock vacuum |0〉 of this field theory is defined by ψ(t, x)|0〉 = 0. Therefore ψ(t, x) and ψ † (t, x) are respectively annihilation and creation operators. The LL Hamiltonian describes the dynamics of N bosonic particles interacting only through a very local potential, a situation which is quite common in cold atom systems. Indeed, the Lieb-Liniger model is a milestone in our understanding of experiments with cold atoms and is nowadays realised in several instances [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18].
• The third and last model is the integrable non-relativistic and classical field theory, called Non-Linear Schrödinger (NLS) model [38,39], with Hamiltonian given by and Poisson brackets {φ(x), φ † ( y)} = ıδ(x − y). This gives rise to the classical equation of motion The NLS equation is a prototypical dispersive nonlinear partial differential equation which, in addition to its application in cold-atom physics, enters the phenomenon of self-focusing and the conditions under which an electromagnetic beam can propagate without spreading in nonlinear media [40]. Properties of the steady state. After a quantum quench, the quantum system undergoes a unitary time evolution. Yet, under quite general assumptions, one expects that the system will locally relax to a certain steady-state which can be characterised in terms of statistical physics. The corresponding ensemble depends though on the nature of the system. If the system under consideration has only its HamiltonianĤ as an extensive integral of motion, it will locally equilibrate to the familiar Gibbs Ensemble (GE) whose temperature is determined by the energy injected at t = 0, which is conserved during the time evolution. On the contrary, if the system is one of those interacting models which are integrable [20,26,37], i.e. those models which, besides the Hamiltonian itself, possess an extensive number of local (and quasi-local) integrals of motionQ j [50][51][52][53][54][55][56][57][58], the late-time steady state will be described by Generalised Gibbs Ensemble [59] (GGE) where the infinitely many generalised temperatures β j are in principle fixed matching the expectation values of the chargesQ j in terms of their value on the initial state. In this case, the presence of infinitely-many conserved charges will deeply constrain the dynamics of these models and many properties of the initial states are not washed out by the time-evolution, but remain encoded in the statistical nature of late-time steady state.

GGE: cahiers de doleances.
It is important to remark that the actual determination of the β j 's from the initial conditions is usually extremely difficult: besides free-to-free quenches (which can be exactly solved by means of Bogoliubov rotations), exact results can be obtained only for very special cases (see Ref. [1] and references therein for a comprehensive overview), leaving out of the analysis most of the possible physically-motivated quenches. From this point of view, any efficient method able to control the GGE, no matter if approximated or based on numerical approaches, would be really welcome. Among the methods used so far, one consists in approximating the GGE by employing only a finite number of charges: this is the idea behind the Truncated GGE (TGGE) [60] and the expectation is that, increasing the number of the charges in this truncation, the corresponding ensemble will converge to the true GGE. While such a method can be successfully implemented in some cases, several pitfalls can undermine its applicability. First of all, one should know the whole set of charges: even though many integrals of motion are often not too hard to be explicitly constructed, determining the whole set of relevant charges is a quest on its own. Indeed, failing to include all the relevant charges led to questioning the applicability of the GGE itself [61], which has been lately redeemed by the discovery of quasi-local charges [50,52,62,63].
Another problem mining the TGGE approach is that the expectation value of the charges could be not-informative at all. A famous example in this context concerns the interacting Bose gas described by the Lieb-Liniger Hamiltonian (2). If one initializes the system in the non-interacting ground state and quenches to the interacting case (κ = 0 → κ = 0), one immediately realises that the expectation values of all the local charges are indeed actually UV divergent [64] (albeit lattice regularizations can be attempted to overcome this issue [65][66][67]), obstructing the applicability of the TGGE. Let's mention that this drawback of the LL model has been lately explicitly overcome through the Quench Action Method [68,69], which relies on the exact computation (in the thermodynamic limit) of the overlaps between the post-quench basis and the pre-quench state [70]. Lately, this method has been used in other models as well (see Ref. [19] and references therein), but its applicability seems to be confined to a restricted set of states [19,71]) while more general (and physical) initial states are still waiting for a solution.
Numerical techniques [72] based on the Quench Action Method have been devised resulting in the ABACUS algorithm [73]: while being able to capture the whole time evolution for quenches with arbitrary interaction changes, its efficiency is drastically reduced for thermal initial states at sufficiently small temperature. In this case, higher temperatures are more difficult to access.
Organization of the paper and main results. In this paper, we are going to revisit the problem through a semiclassical approach, which well approximates the quantum model in the limit of high density and energy, as we thoroughly discuss later on. Semiclassical methods have been applied in a variety of contexts [74][75][76][77][78][79][80] and, for what concerns the Lieb-Liniger model, gave access to predictions unavailable in the purely quantum context. Following the original paper [48], the (semi)classical analysis of integrable field theories out of equilibirum has been extended to weakly-inhomogeneous integrable models [81,82] in the context of the Generalised Hydrodynamics [83,84].
In the following, we focus on the classical counterpart of the Lieb-Liniger model in the repulsive phase (κ > 0), namely the Non-Linear Schrödinger (NLS) equation, and consider quenches from arbitrary homogeneous initial states to the repulsive phase.
• In Sec. 2 we introduce the main character of our analysis, namely the LL model: we briefly review its thermodynamics and its semiclassical limit, i.e. the NLS equation.
• In order to study out-of-equlibrium setups, in Sec. 3 we provide a new analytical toolbox to determine physical observables. More specifically, we determine new close integral expressions for arbitrary moments of the particles' density valid on arbitrary GGEs (see Refs. [46,47] for the quantum result). Furthermore, we also analytically derive exact expressions for the whole probability distribution of the density operator for any GGE, a quantity also known as Full Counting Statistics (FCS). The FCS of several quantities and in a variety of models have been investigated in several instances [85][86][87][88][89][90][91][92][93][94][95][96][97][98][99][100], but for what concerns the LL model only a few results have been obtained so far, despite their relevance in experiments [101][102][103]. For instance, in Ref. [46,47], the FCS for the number of particles in a small interval has been computed in the fully quantum Lieb-Liniger model, at the first order in the size of the interval. On the contrary, in Ref. [75] the FCS for the number of particles for arbitrary intervals has been considered within the semiclassical approximation, but the validity of the result is restricted to thermal ensembles.
• Sec. 4 is devoted to determining the root density from the initial conditions of the out-of-equilibrium protocols. We start presenting the integrability of the NLS from a purely classical perspective and then, building on the findings of Ref. [48], we propose an efficient numerical algorithm for the determination of the GGE. More precisely, the root density is numerically determined evaluating a particular observable, called transfer matrix, which is extensively discussed in this section.
• Sec. 5 is devoted to benchmarking our predictions with first-principle numerical simulations, while in Section 6 we gather our conclusions. A few appendixes follow the main text, providing a few more technical discussions on several aspects touched in this paper.

The Lieb-Liniger model
In this Section, we present some basic facts about the Lieb-Liniger model, in particular we characterise the Hilbert space and the thermodynamics of this model, as well as the expectation values of some physical observables. For the computation of the latter quantities it is worth underlying that the LL model emerges as the non-relativistic limit of the sinh-Gordon model [41,42] or, more generally, of the Toda Field Theories [43]. Given the integrability of the model (see, for instance Ref. [37]), one can look for the common eigenvectors of the whole set of conserved charges. As shown in the original papers by Lieb and Liniger [35,36], these common eigenvectors consist of asymptotically free multi-particle states and their explicit expressions can be found by means of the Coordinate Bethe Ansatz.
Eigenfunctions and eigenvalues. The eigenfunctions can be written as where the rapidities {λ j } N j=1 parametrise the many-body wave-function Above, the sum is over all the permutations P of the rapidities and χ is symmetrically extended when the coordinates are reshuffled in a different order. The coefficients A(P) bear the information about the non-trivial scattering: let Π j, j+1 be the permutation swapping the rapidities at positions j and j + 1, one has with S(λ) the Lieb-Liniger two-body scattering matrix. Eq. (9) makes clear the factorization of the multi-particle scatterings in terms of two-body scattering amplitudes. As we already mentioned, the states (7) are common eigenvectors of all the conserved charges. Besides, the latter act additively on the rapidities, enforcing once again the multi-particle interpretation. For example, for what concerns the energy and the number of particlesN =´dx ψ † ψ one getsĤ where E(λ) = ħ h 2 2m λ 2 is the single-particle energy. In the case of other charges, one gets similar expressions replacing the energy with the proper single-particle eigenvalue of the charge.
Thermodynamics. In order to study the thermodynamic limit of the model, it is crucial to consider a finite system of length L and then take the infinite volume limit, while consequentially rescaling the number of particles in such a way the density n = N /L remains fixed. Hence, we consider the bosons to live on a ring with periodic boundary conditions (PBC): different choices of boundary conditions do not affect the physics in the thermodynamic limit. Analogously to the non-interacting case, where putting the system in a finite volume results in a quantization of the allowed wavevectors, in the interacting case rapidities are constrained by the Bethe-Gaudin equations Apart from the limiting case of free (κ = 0) or hard-core (κ → +∞) bosons, the rapidities are coupled together in a highly non-linear manner and an analytical solution of these equations for large N is hopeless. Fortunately, in the thermodynamic limit we can leave aside the attempt of any exact solution and adopt a coarse-grained approach, which is handled by the Thermodynamic Bethe Ansatz [104] (TBA) method. In the repulsive regime of the LL model, the relevant solutions of the Bethe-Gaudin equations (11) have all real rapidities and, rather than exactly tracking them, we introduce a coarse-grain counting functions ρ q (λ), called the root density. Since later on we will extensively discuss the classical counterpart of the LL model (namely the NLS equation), we find it useful to add a label "q" to denote quantities in the quantum case while, without labels, we will refer to the classical case. Given a physical state as in Eq. (7), Ldλρ q (λ) counts the number of rapidities laying around λ in a small window of length dλ. The root density fully encodes the extensive expectation value of the charges: for what concerns the energy and number of particles, one has and similarly for the other charges. It is however necessary to introduce another quantity, ρ t q (λ), also known as total density, which refers to the total number of available modes in the interval (λ, λ + dλ), while ρ q (λ) counts only the occupied modes which are solutions of (11) and appear in the multi-particle state (7). These two densities are related to each other through the integral equation (11) Together with the root density, it is customary to define the filling and parameterizing it in terms of the effective energy q (λ). In the particular case of a thermal ensemble with inverse temperature β q and chemical potential µ, the effective energy satisfies the non-linear integral equation [104] The importance of the root density. The TBA can be easily generalised to GGEs by including higher charges, provided the whole set of generalised temperatures β j is known [105,106]. Notice that, given a GGE, the root density is fully specified, but the implication holds true also the other way around. Indeed, in Ref. [52] it has been shown that the set of possible GGEs is in one-to-one correspondence with the root densities: in this respect, determining the GGE in the form Eq. (6) is completely equivalent to find the correct root density. This second point of view must be preferred to the previous one, since it does not require to know all the charges explicitly: hereafter, we will always refer to GGEs just specifying the root density (or equivalently the filling). Once the root density is known, one can easily compute the expectation value of the charges (12), but the local relaxation to a GGE implies much more: any local property of the system, such as the expectation values of local observables and correlation functions thereof, is fully specified by the GGE or, equivalently, by the root density.

Moments of the density operator.
While the expectation values of the charges are rather simple, determining other observables (which, in contrast with the charges, undergo a non trivial dynamics after the quench) is in general a hard task. In Refs. [46,47] the problem has been solved for all the moments of the density operator, i.e. 〈(ψ † (x)) n (ψ(x)) n 〉, for arbitrary positive integers n and the expectation value can be taken on an arbitrary GGE (see also Ref. [45,107] for previous results up to n = 4). Since this result will play a key role in the forthcoming discussion of the classical model, it is worth to report it shortly. The expectation values can be obtained expanding in the dummy variable Y the following generating function (we slightly change the notation compared with the original references Refs. [46,47], also in view of the forthcoming semiclassical limits ) where the auxiliary functions ξ n (λ) satisfy the following recursive set of linear integral equa- Above, for an arbitrary test function ω(λ), we define the (quantum) dressing operation as These equations are iteratively solved posing ξ n≤0 = 0 and increasing n step by step. The whole information about the GGE is contained in the filling function ϑ q .

The NLS as the semiclassical limit of the LL
Let us now explore the connection between the Lieb-Liniger model and its classical counterpart, namely the Non Linear Schroedinger equation. We denote with φ(t, x) the classical complex field which obeys the equation of motion (4). As a warm up to understand the applicability of the semiclassical approximation, it is useful to look at thermal ensembles of the LL model without referring to its integrability structure. As we have already commented in the introduction, the semiclassical limit can be viewed as a high temperature limit. In order to see this, we define a rescaled Hamiltonian Comparing with the Hamiltonian (2), we haveĤ = ħ h 2 2mĤ , hence the equality of the density matrices with As it is clear, the ħ h → 0 can be interpreted as a high temperature/small coupling limit of for the rescaled Hamiltonian: we will pursue this second interpretation along this section.
Partition functions. Let us initially consider the (quantum) partition function for the LL model. We consider the rescaled Hamiltonian (20) with an inverse temperature β q (22) and a chemical potential µ, describing the partition function within a path integral approach.
There is a standard procedure to reach the semiclassical limit of this expression (see, for instance [108]) and this consists of supplementing the real space with a Euclidean time coordinate τ. Assuming periodic boundary conditions in the real space, the quantum problem is mapped into a classical partition function on a torus (x, τ) Fig. 2): for historical reasons, we refer to this domain as the Matsubara torus. The quantum partition function can then be written as where ψ(τ, x) now denotes a classical field whose domain is the Matsubara torus. On the other hand, the classical partition function can be written as a path integral on a ring x ∈ [0, L] (see Fig. 2) In order to map the quantum partition function into a classical one, we should be able to "shrink" the imaginary time dimension β q → 0 (see again Fig. 2): this is indeed what the ħ h → 0 limit does, as it is clear from Eq. (22). To do so, let us consider the Fourier components of the field along the Euclidean time direction Above, for later use, we have suitably normalised the modes: this will introduce an overall constant renormalization of the partition function which however does not affect the thermodynamics. This is dominated by the zero-frequency mode and therefore it is described by the classical NLS. Indeed, plugging this mode decomposition in the quantum action, one finds and, when β q → 0, this path integral is dominated by the zero-frequency Matsubara mode since the (i2π/β q )nψ † n ψ n term tends to pin ψ n =0 → 0, while it vanishes for ψ 0 that is then free to fluctuate. Hence, discarding the non-vanishing Matsubara frequencies we can approximate the partition function as which is nothing else than the classical partition function Eq. (24), provided the replacement Notice that the classical inverse temperature β is kept constant in the ħ h → 0 limit, as it should be.
TBA in the semiclassical regime. After this warm up on the partition functions, let's now move on to discuss the semiclassical limit from the integrability point of view within the TBA framework, extending the original discussion done in Ref. [48] (see also Ref. [83]). We will now use the quantum TBA to access the classical one, together with other results. Once the mapping between the quantum and the classical TBA has been established, in the following we will refer to the classical TBA simply as the TBA of the NLS model. As already stated, the classical physics emerges from the quantum one in the limit of high occupation numbers, which in free systems is the mode density: within the interacting integrable framework, the mode density is replaced by the root density, which is then expected to diverge with ħ h. Inspired by the thermal case and using Eq. (25), one immediately sees that in the semiclassical limit the density of particles diverges as ∝ (2mκħ h 2 ) −1 : this same behavior is reproduced at the level of TBA if one poses with ρ(λ), interpreted as the classical root density, kept constant while sending ħ h → 0. Let us use this scaling to extract the semiclassical limit of the basic building blocks of the TBA: we will be back to the thermal state later on. We start defining the classical total root density: in this respect we consider Eq. (13) with the replacement (29). As it stands, the left and right hand sides do not have a well defined ħ h → 0 limit: using the identity ρ(λ) =´dλ δ(λ − λ )ρ(λ ) we can recast Eq. (13) in the following equivalent form . By means of a simple integration by parts, we obtain the following equivalent form, which is better suited for taking the semiclassical limit Taking the ħ h → 0 limit one has where the singular integral is handled with the principal part prescription. We can now define the classical total root density trough the limit which defines the total root density in the classical model. In analogy with the quantum case, one defines the classical filling ϑ(λ) = ρ(λ)/ρ t (λ): this definition, together with Eq. (33), is consistent with the one already given in Ref. [48]. A similar analysis can be performed on the dressing (19) and on the equation for the thermal state (15): the classical dressing on a test function ω(λ) is defined as and one has the following quantum-classical correspondence Concerning the thermal states, one can still obtain convenient integral equations in terms of an effective energy (λ), although its relation with the filling is modified with respect to the quantum case (14) Then, the thermal state is determined by the following integral equation which can be readily derived from Eq. (15) in the ħ h → 0 limit, as done in Eq. (28). The comparison between the expression of the classical filling Eq. (36) and the quantum one Eq. (14) in terms of the respective effective energies deserves further comments. Indeed, Eq. (36) is the generalization to the interacting integrable case of the Rayleigh-Jeans distribution, which is obeyed by the classical fields at equilibrium and leads to the famous UV catastrophe, namely to an UV divergence of the energy expectation value 〈H〉 =´dλ E(λ)ρ(λ). Indeed, for large rapidities, Eq. (37) gives (λ) ∼ β E(λ) and, in view of Eq. (36), one readily gets ρ(λ) ∼ (β E(λ)) −1 with a consequent UV divergence of the energy. The UV catastrophe is what led Planck to discover quantum mechanics: in the quantum case Eq. (14) the Fermi-Dirac distribution gives a finite energy and regularizes the classical prediction. The UV divergence of the local charges on a physical state, such as a thermal ensemble, suggests that they are not the best observables to look at. Furthermore, being our main interest out-of-equilibrium protocols, we wish to study observables which have a non-trivial time evolution after the quench: in this way, we probe the time evolution and observe the relaxation to the GGE. In view of all these considerations, we focus our attention on the moments of the density |φ(x)| 2n : we now derive close integral expressions valid on arbitrary GGEs from the semiclassical limit of the analogue expression in the LL model [46,47]. This is a new result in the context of the Non-Linear Schrödinger equation.
The density moments from the semiclassical limit. The semiclassical limit of the generating function Eq. (16) is straightforward. Indeed, from the analysis of the thermal case we know which immediately leads to the classical generating function (39) Above, we plugged the correspondence Eq. (35) in the definition of G q n (16) and defined As the last step, we plug the above definition of ζ n in the integral equations for ξ n Eqs. (17) and (18) and take the semiclassical limit, resulting in the desired equations for the classical case Within the quantum case, the expression for the density moments of Ref. [46,47] has been found to be in agreement with previous results known in the literature up to the fourth moment [45,107]. However, continuous quantum models are notoriously difficult to be simulated and these results lacked any numerical benchmark: the situation is different in the classical case, where the thermal expectation values can be efficiently computed by means of the Transfer Matrix (TM) method. We leave to Appendix A the description of the TM approach, while in Fig.  3 we benchmark the TBA results for different choices of temperature and chemical potential. Lately, in Section 5, we will use the density moments as a diagnostic tool to compare GGE's The Full Counting Statistics. Apart from the moments of the density 〈|φ(x)| 2n 〉, one could ask what is the probability to measure a certain value of |φ(x)| 2 . This is especially important when one can access a limited number of samplings or it is interested in the extreme-value statistics. Let us introduce P(d) as the probability, on a given GGE, that a measurement gives |φ(x)| 2 = d as an output: formally, we can express it as the average of a Dirac−δ as follows Above, the expectation value is taken with respect to a given GGE and the Fourier representation is introduced for later use. In the literature, the probability P(d) is often called Full Counting Statistics (FCS) (in this case, FCS of the density operator) and it is notoriously hard to be analytically accessed. Here, we present an exact result for the classical FCS of the density operator valid for arbitrary GGEs. Therefore, at the price of leaving aside finite intervals, we can generalise the result of Ref. [75] to out-of-equilibrium situations. In this case, an analogue result in the quantum case from which we can extract the semiclassical limit is not available and one must proceed from scratch. The starting point is an expression similar to Eq. (39), but of wider generality ∞ n=0 q 2n (16) n (n!) 2 〈|φ| 2n 〉 = exp where In principle, one could recover Eq. (39) power-expanding in q the r.h.s. of Eq. (43) (and the integral equation for W q as well). However, the other way around is far from being trivial: we derive Eq. (43) from known results in the classical sinh-Gordon model and taking the non relativistic limit. This derivation is rather technical and it is provided in Appendix B, while here we build on Eq. (43) in order to access the FCS. The analysis presented in Appendix B resulting in Eq. (43) assumes real parameters q, however in order to access the FCS we need to perform a proper analytical continuation. Let us introduce an auxiliary function F (p) which has the following power expansion We note that this is the same power expansion of Eq. (43) if one replaces q 2 → p. Using the function F (p), we can reach a convenient expression for the FCS. Indeed, from Eq. (45) one readily has and establishes the following chain of equalities Above, we regularize the expression with an infinitesimal shift in the complex plane γ → γ+i0 + . Using this result in the Fourier representation of the full counting statistics (42) and performing the integrals in k and γ, one reaches the following representation of the FCS Above, J 0 (x) is a modified Bessel function of the first kind. As it is clear from the above, the knowledge of the F (p) function gives access to the FCS by means of a simple convolution. Going back to Eq. (43), one sees that the expression as it stands gives access to F (p > 0) (as it is readily seen from Eq. (45)). However, extracting the FCS requires the knowledge of F (p < 0): we fill this gap performing a proper analytical continuation of the r.h.s. of Eq. (45). While doing so, some subtleties arise. Let us proceed naively and promote q to be imaginary in Eq. (43) posing q = iτ, furthermore we define the auxiliary function s τ (λ) = ϑ(λ)W iτ (λ): by means of a direct replacement in Eq. (45) one finds where s τ (λ) satisfies From the left hand side of Eq. (49), we see that the same result on the r.h.s. must be obtained if one replaces τ → −τ. By a direct inspection of Eq. (50) we readily see s −τ (λ) = s * τ (λ), therefore after the replacement τ → −τ in Eq. (49) we get Imposing the consistency of the above with Eq. (49), we can conclude that F (−τ 2 ) must be real (which was expected, since the FCS must be a real function) and we must necessarily havê This equality can be easily shown as follows. We use a vector-matrix notation introducing a basis of states |λ〉 and defining an operator Ω(τ) with matrix elements We then introduce the state |s τ 〉 such that s τ (λ) = 〈λ|s τ 〉 and analogously the state |1〉 such that 〈λ|1〉 = 1: matrix-vector contractions and scalar products are computed integrating over the |λ〉 states. In this notation Eq. (50) can be written as Above, we used the hermicity of Ω(τ) and introduced a complete basis of eigenvectors Ω(τ)|n τ 〉 = µ n (τ)|n, τ〉. In the vector-matrix notation we havê which is clearly real, since Ω(τ) has real eigenvalues. This proves Eq. (52). There is another subtlety in the analytic continuation, as it can be seen from Eq. (55): if one of the eigenvalues of the operator Ω(τ) vanishes, the overlap 〈1|s τ 〉 develops a first order pole singularity as a function of τ. Strictly speaking, such a singularity is not integrable and causes Eq. (49) to be ill-defined. This singularity arises in the analytical continuation while q = iτ approaches the imaginary axis: if we move slightly off from the imaginary axis the operator Ω(τ) (53) is no longer real and, for infinitesimal shifts, the eigenvalues µ n acquire a small imaginary part. This suggests to regularize Eq. (55) replacing µ n (τ) → µ n (τ) ± i0 + , both of sign choices leading to the same result, as we discuss below. Indeed, let us use the (regularized) expression (55) into Eq. (49), which results in Figure 4: Comparison between probability distribution for the particle number density P(d) ≡ (|φ| 2 = d) on thermal states computed extracting the root density from TBA equations (37) and using the Transfer Matrix approach. Here β is the inverse temperature characterising the state.
Above, we summed and substracted the singular part, letting τ j be the zeroes of the n th j eigenvalue µ n j (τ j ) = 0. The first term is regular and real, the second term is easily integrated leading to (we focus on the contribution of a single pole) . (57) The above expression is not explicitly real and has a non-trivial dependence on the sign ±. Both issues are resolved if it holds true for any τ j , then which is real and analytic and does not depend on the sign of the regularization. The identity Eq. (58) can be analytically proven, as we show in Appendix C. Employing Eq. (58), we are led to a well-defined expression for the F function (60) In Fig. 4 we evaluate the exact result (48) on thermal states and compare it with the numerical data obtained with the TM approach, finding perfect agreement. However, let us stress once again that Eq. (48) is of much wider applicability and holds true on arbitrary GGEs, hence can be used in out-of-equilibrium protocols as well, as we will do in Sec 5.
The semiclassical limit in the lab. Along this paper, we resided to the small parameter ħ h as a convenient way to take the semiclassical limit. While this is probably the most natural approach, on the other hand it is not of immediate experimental applicability: however, as we have already commented at the beginning of this section, the same limit can be analogously obtained in a high temperature and weak coupling limit with ħ h kept constant, where the quantum temperature and the interaction are sent to zero while keeping their ratio fixed. More specifically, let us consider the Lieb-Liniger Hamiltonian Eq. (2) with the following parametrization and let us consider a thermal state e −β q (Ĥ−µN ) , then the classical thermal state Eq. (24) is achieved in the β q , g → 0 limit, provided we keep the ratio fixed in such a way that the classical temperature β is fixed β = β q /((2mħ h −2 )g). In this case, the quantum expectation values go to the classical ones through the dictionary ψ → φ/ 2mħ h −2 g. The semiclassical limit of thermal states and the corrections induced by quantum fluctuations are well controlled (see eg. Ref. [109]) and out-of-equilibrium setups in the weakly interacting and high density regime are expected to be well described by the classical physics as well.

The classical integrability of the NLS: the inverse scattering method
In this section, following Refs. [38,39,49], we provide a brief summary of the classical integrability of the NLS equation. The main aim is to put in evidence the connection with the semiclassical limit of the quantum problem and to find an efficient procedure to extract the classical root density from the initial state. For certain initial conditions with zero energy density, see e.g. Refs [110,111], the inverse scattering problem can be completely solved and the whole time evolution predicted. On the other hand, as we extensively commented, our main focus is on initial conditions with a well-defined thermodynamic limit, i.e. with a finite energy density, and in the determination of the final steady state.
Zero curvature condition. The remarkable observation which led to the development of the inverse scattering trasform is that the evolution equation can be seen as a compatibility condition for a linear scattering problem depending parametrically on the spectral parameter λ. More specifically, for the NLS model, one introduces the following system where U λ , V λ depend on the spectral parameter λ and on the field φ(t, x) = u(t, x) + ıv(t, x). Expressed in terms of standard Pauli matrices, they read (we omit the explicit space-time dependence to simplify the notation) The two-component field F (t, x) can be expressed in terms of the eigenstates of σ z as Notice that the system of the two differential equations (62a) and (62b) are compatible only if it holds the so called zero curvature condition This equation has a clear geometrical interpretation in terms of gauge fields. Indeed, let's define the covariant derivative as where Consider the differential form A = A µ dx µ = V λ dt +U λ dx and a closed path Γ in space-time. In light of Eq. (65), the differential form A is a closed form and an application of Stokes theorem gives P exp where P is the path ordering operator. The above operator implements the parallel transport, which is trivial on a closed path because of the vanishing curvature.
Monodromy matrices. Let's now consider the closed path defined by the points (t 1 , x), (t 1 , x), (t 2 , x), (t 2 , x) and (t 1 , x) (see Fig. 5) and define the propagators (monodromy matrices), Using (68) and (69) we can write, From the elementary properties S −1 λ (t 1 , t 2 ; x) = S λ (t 2 , t 1 ; x) and T −1 Hence, if there exist two points in the space-time such that then, in view of the relation (71), the quantity is time-independent. The simplest ways to satisfy Eq. (72) are to take Note that, for a given configuration of the field φ(t 0 , x) at given time t = t 0 , the propagator T λ (t 0 ; x, y) can be efficiently computed numerically by simply discretizing the path-ordered integral in (69) and expressing it in terms of the following matrix product Expanding Eq. (71) for t 1 → t 2 , one sees that T λ (t; x, y) and M (λ) = V λ (t, x) form the socalled Lax pairṪ Then for any value of the spectral parameter λ, the quantity defined in Eq. (73) can be used to generate an infinite set of conserved quantities.
Thermodynamic limit. As explained in Sec. 2, in the thermodynamic limit the stationary states of the Lieb-Liniger (and correspondingly of the NLS) model is completely described in terms of the root density. We already mentioned that via Eq. (12), the root density is in oneto-one correspondence with the expectation values of the conserved quantities. Since, in the semiclassical limit, the conserved quantites are all generated by Eq. (73), it is natural to expect a relation between T λ (x, y) and the root density. In order to unveil it, we first consider the system on a finite volume of size L and then take the limit L → ∞ at finite energy density. We thus enforce periodic boundary conditions for x ∈ [0, L]: then Eq. (71) implies T λ (0, L) is conserved for any λ. By performing a series expansion at large λ [38,48,112], one can show indeed that where Ω(λ) is the generating function of all local conserved densities via Now, since σ x U λ σ x = U * λ , the diagonal entries of T λ (0, L) are complex conjugate. Combining this fact with Eq.  Figure 6: Root density generated in two different quenches (continue lines) compared with the thermal state it would be reached in the absence of integrability (dashed lines). More specifically, the system is prepared in a non-interacting thermal state with Hamiltonian H t<0 =´dx |∂ x φ| 2 and inverse temperature β = 1 and densities chosen to be in the strong-interacting regime D = 〈|φ| 2 〉 equal to 1 and 2 respectively. The parameters of the putative thermal states shall be chosen to match the initial value of energy and density. Note that because of the UV divergence of the energy, the final value of β = 1 remains unchanged, while µ is adjusted to match the initial densities. Inset: tails of the roots in logarithmic scale, which besides of collapsing on the same curve they also agree with the initial thermal distribution, as we explain in the main text.
Post-quench root density. In the quench protocol we have considered in this paper (see next Section for more details), the system is initially prepared in a thermal state and we are interested in the post-quench root density which provides its GGE representation. This can be achieved 1. drawing a field configuration φ(x) from the Boltzmann measure given in Eq. (24); 2. computing the matrix T λ (0, L) via eq. (74); 3. using Eq. (76) to extract the root density as Note that the root density is self-averaging and therefore for sufficiently large L, a single configuration would be sufficient. In practice, for the numerical implementation, we have averaged Eq. (79) over many field configurations so that finite size fluctuations were suppressed (see below).

Benchmark of the analytic results: time evolution vs GGE predictions
In this section we compare the analytical findings of the previous section with ab-initio numerical simulations of the NLS (details in App. D). In the quantum quench protocol, the initial state is often chosen in the form of a pure state in order to enhance its quantum properties, but obviously density matrixes (such as thermal states or GGEs) are equally physically-motivated states. In the perspective of using the classical theory as an approximation of the quantum model, we have initialized the system directly in a Gibbs ensemble of the NLS and we have let the system evolve with a different interaction for t > 0. In order for the GGE prediction to emerge, an averaging procedure is required: within the quantum context, the expectation value on the initial state accounts for both quantum and statistical fluctuations, while in the classical case the first are obviously absent. Hence, we average with respect to the initial random configurations. More specifically, the microscopic protocol consists in three steps: • Sample the initial field configuration from the desired thermal ensemble.
• Evolve each field configuration with the deterministic NLS equation, measuring the observables along the evolution.
• Repeat the above steps for several initial field configurations and average over them.
As stated before, in the thermodynamic limit, the spatial average of a single field configuration would be sufficient. In practice, as we work at finite volume, we average over the initial ensemble to further suppress the fluctuations. The numerical computation of the root density on the initial state, which is then fed to the expressions of Section 3 to access the steady-state's observables, is performed in similar steps. For any field configuration dragged from the initial ensemble, we compute the associated transfer matrix through a brute-force computation of Eq. (74), then the root density is extracted by means of Eq. (79) computed at finite volume: in order to extract the thermodynamic limit, we use the self-averaging property of the system and average the root density over the initial configurations (see App.D for further details). The averaged root density is then used in computing the moments of the density and the whole FCS. In the following, we start from non-interacting thermal states based on H t<0 =´dx |∂ x φ| 2 and then evolve with Eq. (4). Thermal states in free systems (and, more generally, GGEs) are easily generated in the Fourier space, where the modes are independently gaussianly distributed (see App. D). The sampling of an interacting thermal states can be achieved with a Metropolis-Hasting algorithm [113], similar to what it has been done in Ref. [81]: however, we do not expect any qualitatively-different behavior, hence we consider the non-interacting case for the sake of simplicity.
In Fig. (6) we compare the root densities obtained through different quenches. More specifically, the system is initialized in a non-interacting thermal state and then let evolve in the presence of interactions. We point out that at large rapidities the root density is expected to coincide both with the putative interacting thermal state, as well as the non interacting thermal state. This can be forecast in view of the fact that large rapidities are weakly interacting, hence they are not affected by changes of the interaction.
In Fig. 7 we show the time evolution of the density's moments for different choices of the initial state, finding perfect agreement with relaxation to the forecast GGEs. Instead, in Fig.  8 we focus on the time-evolution of the full counting statistics: the time-evolved distribution clearly approaches the GGE prediction, confirming once again the validity of our method. The non-monotonous behavior of the FCS unveils the highly interacting nature of the system. Indeed, in the absence of the interactions and thus on a gaussian steady state, the FCS is straightforwardly shown to be a simple exponential P(d) = D −1 e −d/D (with D the average density): in Fig. 8 at time t = 0 the system is initialized in a gaussian ensemble and the curves are perfect exponentials. Then, at later times the effect of interactions becomes more and more important, shifting the peak of the distribution from d = 0 to a finite value. At equilibrium and for finite intervals, this features have been extensively studied, for example, in Ref. [75].
This classical out-of-equilibrium protocol can be viewed as an approximation of a quench in the quantum Lieb-Liniger, from zero to a small value of the interaction and starting from a high temperature thermal state. Quenches starting from weakly interacting initial thermals states, which are approximated by interacting classical thermal ensembles, can be considered as well and are a simple application of our general method.  Fig.6 and we focus on the D = 2 case. Since the initial state is gaussian, at time t = 0 we have ω K = 1, however as time passes the strongly interacting nature of the system is unveiled and the values of ω K depart from unity, relaxing to the predicted asymptotic value with excellent accuracy. In our microscopic simulation we experienced important finite-size effects and the data shown in the plot are extrapolated to the thermodynamic limit. (b) We provide a brief analysis of the finite-size effects focusing on the time evolution of the first non trivial observable ω 2 . Inset: extrapolation to the thermodynamic limit of the plateau extracted from the microscopic evolution (blue line) and of the plateau obtained from the root density (yellow line). In both cases, we observe a clear ∝ L −1 approach to the asymptotic values, which are in perfect agreement.

Conclusions
In this paper we thoroughly analyzed quench protocols in the classical Non-Linear Schröedinger equation, which can be viewed as the semiclassical limit of the interacting Bose gas. Our interest was twofold: on the one hand, we aimed to study the classical model taking advantage of existing results on the quantum model, projected in the classical world through a proper semiclassical limit. To this end, we determined new exact expressions for the moments of the local density in the form of recursive integral equations, which are valid for arbitrary GGEs. These expressions are the classical analogue of the study of Ref. [46,47]. Moreover, viewing the classical NLS as the non relativistic limit of the classical sinh-Gordon model, we provided the exact probability distribution of the local density: this new result has no quantum analogue so far.
On the other hand, we were interested in investigating the quantum model within the semiclassical regime: a key-point of addressing a quench protocol is being able to determine the emergent GGE from the initial conditions and, within the quantum case, there are no analytical or numerical methods to answer this question (with a few remarkable exceptions [70,72]). We completely solve this problem in the classical regime, providing an efficient numerical procedure to determine the GGE from the initial state.
Several interesting questions naturally emerge from our investigation. First of all, even though we mainly focused on the NLS, the same approach can be applied to study other classical integrable systems. In particular, it would be interesting to address systems whose TBA possess more than a single species of excitation (string) and lattice systems, such as the Toda lattice [38]. It would be interesting to gain further insight in the semiclassical limit of other relevant quantum models, such as the XXZ spin chain and its classical counterpart, the Landau-Lifshitz model [114].
Hybrid semiclassial-quantum methods could be imagined as well, with the semiclassical computation of the initial root density, which can then be used as an initial condition for quantum methods, for example the generalised hydrodynamics [83,84]. This could be implemented, for instance, to study the the famous Newton-Cradle experiment [9]: the state after the initial Bragg pulse still lacks a faithful description [115].
Lastly, it would be very intriguing to improve the semiclassical approximation of the quantum model: the classical system can be regarded as the zeroth order in the ħ h expansion of the quantum one. Is it possible to include further orders in ħ h, providing an expansion of the quantum model building on the classical theory? Can we include these quantum corrections in our method to determine the root density? We leave these exciting questions open for future investigations.

A The transfer matrix method
The transfer matrix method takes advantage of the classical/quantum correspondence of the partition functions to map the 1d classical problem in a 0d quantum one, reducing the problem of computing classical thermal averages of local observables to find the ground state of a two-dimensional Hamiltonian [48,[116][117][118][119]. For the sake of clarity, we consider directly the momenta of the density, but the same procedure can be repeated for any observable. The classical expectation value is given by Lately, we will consider the limit L → ∞. We rewrite the above as a quantum expectation value introducing a two-dimensional quantum-mechanical Hilbert space H T M spanned by wavefunctions having support on (u, v) ∈ 2 , together with the quantum Hamiltonian H T M Then it is a straightforward excercise to establish the correspondence between the classical expectation values and those in the ancillary quantum problem where | j〉 are the eigenvectors of the quantum Hamiltonian H T M | j〉 = E j | j〉. Above, the trace appears because of the spatial periodicity of the classical model: the length L plays the role of an effective inverse temperature for the auxiliary quantum problem. If we are now interested in the thermodynamic limit we let L → ∞, hence the quantum trace is projected on the ground state. Thus, we finally have where Ψ G is the ground state wavefunction of the Hamiltonian H T M . Rather than approaching directly the two dimensional problem, it is convenient to use the rotation symmetry in terms of polar coordinates u = r cos θ and v = r sin θ The ground-state wavefunction is symmetric under rotations, hence it is independent on θ and as such it can be found (numerically) solving the problem in the single-variable r.
B Derivation of Eq. (43) In this appendix we derive Eq. (43), which is the first building block to access the FCS. The starting point is the classical ShG model: in contrast with the quantum Hamiltonian (1), we remove the ħ h coupling and set the mass scale m 0 = 1/2. Therefore, we consider the classical Hamiltonian where Π and Φ are classical conjugated fields {Φ(x), Π( y)} = δ(x − y). Within the quantum case, the non-relativistic limit has been thoroughly examined to study the LL model starting from known results in the ShG field theory [41,42,120]. Indeed, the result for the onepoint functions of Ref. [46,47] has been obtained taking the NR limit of a result of Negro and Smirnov [121,122] (see also Ref. [123]) in the quantum ShG. Here, we follow a similar procedure, albeit in the classical case.
The classical analogue of the Negro-Smirnov formula has been derived in Ref. [81] and it consists in a recursive set of integral equations for the expectation values of certain operators in the ShG model, the formula being valid for arbitrary GGEs Above, ϑ ShG (θ ) is the filling fraction in the classical ShG theory and p k (θ ) satisfies Presenting a detailed analysis of the NR limit goes beyond the aim of this work, which is mostly focused on the semiclassical case: we leave to the original references [41,42] the details (even though they concern the limit of the quantum model, the analysis in the classical case is exactly the same) and quote the key relations. The NR limit is attained sending c → ∞, while instead g → 0. In particular, the NLS model (3) is attained if one sets At the level of TBA, one connects the relativistic rapidity θ and the Galilean one λ looking at the momentum eigenvalue λ = c 2 sinh θ : imposing λ to be finite in the c → ∞ limit, one must consider small θ −rapidities, therefore at first order one simply has θ 2λc −1 . Finally, the filling in the ShG simply goes to the NLS one, provided we correctly link the rapidities ϑ(λ) ϑ ShG (2λc −1 ) .
We can now proceed and take the NR limit of Eq. (86), from which Eq. (43) will follow. Rather than naively take c → ∞ (which would lead to a trivial result), we first rescale q = kc −1 and take the NR limit keeping q fixed. Looking at the l.h.s. of Eq. (86) (and rescaling the coupling g as per Eq. (88)), one gets We retain up to first order in c: the O(c 0 ) term in the above exactly cancels the first term in the r.h.s. of Eq. (86) and one should compare the next order. By means of a direct power expansion and using Eq. (90), one readily gets lim c→∞ 〈e q4Φ 〉 = n q 2n (16) n (n!) 2 〈|φ| 2n 〉 , which is the l.h.s. of Eq. (43). The r.h.s. is readily found taking the NR limit of Eq. (86) (and Eq. (87): while doing so one defines W q (λ) = lim c→∞ (p cq (2λc −1 )), then integrating in q and inverting the logarithm appearing in Eq. (91). Eq. (43) then readily follows.

D.1 Solution of the equation of motion
We now consider the problem of solving the equation of motion starting from a given initial field configuration We approach the problem through a Hamiltonian-splitting procedure [124], which is stable and exactly conserves the number of particles. We split the Hamiltonian H in three parts with The field is discretized as φ j,i ≡ φ( j∆t, i∆x) with i = 0, . . . , N − 1 and j = 0, . . . , M − 1. Then, each step in the time evolution is performed through three different phases The first and third steps are easily performed in the real space, while the one in the middle is performed in the Fourier space, where H 2 acts diagonally. φ j+2/3,k = e −i∆t E kφ j+1/3,k , E k = 2 ∆x 2 1 − cos In the limit N → +∞, ∆x → 0 with L ≡ N ∆x fixed, defining λ ≡ 2π L k, we have E k → E(λ) = λ 2 at leading order.
The energy scale, dominated by the kinetic term, sets the allowed values of ∆t, which must always fulfill ∆t max In our simulations, we always used (∆x) −2 ∆t 10 −3 . In order to extract the continuum limit, we extrapolate the finite-size data according with a simple Taylor-expansion ansatz in the length of the system

D.2 Generation of initial states and averages
As explained in the main text, for the sake of simplicity we focus on quenches from the free theory to the interacting one. Free thermal states (and, more in general, GGEs) are simply formulated in the Fourier space, where the modes are independent and gaussianly distributed with zero mean. Hence, field configurations are easily generated in the Fourier space and then converted back in the real one. On a thermal state with inverse temperature β and chemical potential µ, the modes are distributed according with the Raylegh-Jeans distribution The density of particles D fixes the chemical potential according to the relation (valid in the continuum limit) Given an initial configuration C the value of a generic observable O depends on it. The ensemble average is estimated as and analogously for the standard deviation.

D.3 Extraction of the root density
The key identities for the extraction of the root density are Eqs. (74)- (79), which explicitly read where τ = 1 2 |φ| 2 − λ 2 . For a system of length L, the matrix T λ is approximated as, W ∆x (t, i∆x), N ∆x ≡ L, N 1 (112) and the root density is easily extracted. The root density is then used as an input for the exact formulae for the one-point functions Eq. (41) and the FCS Eq. (50). A case amenable of a fully analytical check is when the initial field configuration is homogeneous in space. In this situation the matrix T λ does not depend on x and a straightforward calculation shows that, which is a semicircle distribution, as we show in Fig. 9 (a). In Fig. 9 (b) we consider the numerically-extracted root density choosing as initial states free thermal states of different temperatures. We notice that this simple result is in agreement with the analytical solution of the quantum quench from the non-interacting ground state. Indeed, the root density computed in Ref. [70] collapses to the semi-circle low in the limit of small interaction, i.e. where the classical model emerges.