Generalized hydrodynamics of the attractive non-linear Schrӧdinger equation

The relentless experimental advances in probing strongly correlated quantum matter are hardly overemphasized [1, 2] and triggered an as much impressive amount of theoretical research. The one-dimensional world stands out as a fruitful melting pot for experiments [3–7], numerical methods and analytical expressions: indeed, particles moving on a line cannot avoid each other, with the consequent enhancement of interactions and correlations. As a consequence, the fascinating problem of describing nonequilibrium phases of matter in low dimensions is extremely challenging and often far beyond the comfort zone of perturbative methods. Hence, other approaches should be envisaged.

Tensor networks [8] are one of the most versatile and powerful numerical methods to access one-dimensional physics, but they come with their own limitations: the nonequilibrium late-time behavior of many-body system is extremely hard to capture, due to its highly entangled nature. If on top of this one also wishes to set aside lattice systems and explore continuous setups, the problem becomes quickly intractable. Other strategies should then be used, for example semiclassical methods or exact solutions: this paper lies at the interface between these two.

A widely known modelling of 1D cold atom experiments is through the interacting Bose gas Hamiltonian, which in dimensionless units reads

$$\begin{equation}\hat{H}=\int \mathrm{d}x\enspace {\partial }_{x}{\hat{\psi }}^{{\dagger}}{\partial }_{x}\hat{\psi }+c{\hat{\psi }}^{{\dagger}}{\hat{\psi }}^{{\dagger}}\hat{\psi }\hat{\psi }+V(x){\hat{\psi }}^{{\dagger}}\hat{\psi },\end{equation} \tag{ 1 }$$

with c describing a local pairwise interaction and V(x) being the shallow confining potential and $[\hat{\psi }(x),{\hat{\psi }}^{{\dagger}}(y)]=\delta (x-y)$. This model hardly needs any introduction due to its widespread use in both experimental [9–19] and theoretical communities.

In experiments, interacting bosons squeezed in elongated traps, either by means of optical lattices or atom chips, are trustfully described by this Hamiltonian. Remarkably, both the trapping potential and the interactions are highly tunable parameters, also in real time. Particularly interesting for us, as we discuss later on, is the possibility of manipulating the interaction by means of Feshbach resonances [20–23].

The ground state and low-energy properties of the repulsive (c > 0) 1D Bose gas can be efficiently described within bosonization [24, 25], but depending on the question at hand, the usefulness of this approach is at stake: macroscopic changes of the trap or interactions will create excitations high in the spectrum, far beyond bosonization applicability. In the opposite limit, semiclassical methods take the lead: in the case of large mode occupancy, one can discard the quantum nature of the fields by replacing them with classical variables $\hat{\psi }(x)\to \phi (x)$. The system is now described by the classical non-linear Schrӧdinger (NLS) equation [26, 27]

$$\begin{equation}\mathrm{i}{\partial }_{t}\phi =-{\partial }_{x}^{2}\phi +2\varkappa \vert \phi {\vert }^{2}\phi +V(x)\phi (x),\end{equation} \tag{ 2 }$$

where ϕ is a classical complex field and ϰ the classical interaction strength. This approximation holds when interactions are weak and curiously emerges in two opposite scenarios. On the one hand, it can describe a single macroscopically occupied mode (quasicondensate) and within this setup is often referred as Gross–Pitaevskii equation [28]. Since this description captures a quasicondensate, it can be interpreted as a low temperature limit.

On the other hand, high temperature limits of the Bose gas can be treated with the NLS as well, which now describes a fluctuating classical field [29–31]. The fact that cold atoms can be described by means of a high temperature fluctuating classical field could appear as a contradiction, but the concept of 'high temperature' must be compared with the energy scales of the 1D experiments, i.e. the interactions and the oscillator frequency of the longitudinal trap, which can be rather weak. The appeal of semiclassical approximations is self-evident: classical fluctuating fields can be easily simulated with Monte Carlo methods on very large scales and for long times. Besides, external perturbations and noise are easy to implement.

Outside the semiclassical realm, whose applicability is not restricted to the Bose gas, one can probe the true quantum nature of the system with analytical methods. Indeed, the 1D Bose gas (in the absence of the trap) is an outstanding example of an integrable model, originally solved by Lieb and Liniger in 1963 [32, 33] and often referred to as Lieb–Liniger (LL) model after them. When applicable, integrability is powerful: it gives access to exact results directly in the thermodynamic limit, ranging from equilibrium physics [34] to out of equilibrium setups [35]. As a downside, deriving analytical solutions can be challenging and until recently the effects of soft integrability breaking terms, such as the longitudinal trap, were hard to include in the theory.

The table turned in 2016 with the proposal of a new hydrodynamic approach to integrability [36, 37] now dubbed as generalized hydrodynamics (GHD), which is built only on the local and instantaneous integrability of the Hamiltonian, making it the ideal tool to describe several experimental setups. The interested reader can refer to the recent reviews [38–44]. Shortly after the original papers, inhomogeneities in the Hamiltonian such as external potentials [45] and time-changes of the interactions [46] were readily included. Recently, GHD predictions have been challenged by real experiments showing a spectacular agreement between the two [40, 47–49]. Noticeably, the very defining properties of integrability is having an extensive number of local conserved quantities that hinder thermalization and cause relaxation to generalized Gibbs ensembles (GGE) [50], as probed in early experiments [51, 52], and allows for novel nonequilibrium phases. Very interestingly, integrability can stabilize quantum matter in regimes that simply do not exist in its absence: the Bose gas within the attractive phase c < 0 is an outstanding example of it. At equilibrium, the attractive Bose gas is extremely unstable [53–55]: in this phase, the gas can form bound states [56–59] and for any finite temperature the gas collapses in a gigantic molecule consisting of all particles. However, integrability overturns the rules of the game by stabilizing the few-particles bound states and allowing stable nonequilibrium phases. This has been already realized and discussed in experiments [22, 23], but theoretical results are scarce [60, 61] and for a good reason: one should initialize the system in a far from equilibrium setting and this is extremely difficult to be analytically addressed.

Fortunately, GHD comes in handy when facing this challenge: in a recent work [62], the three of us proposed to initialize the gas in the attractive regime by means of a slow interaction tuning starting from the repulsive phase. Within each phase, GHD describes the interaction changes [46] and must be supplemented by a description of what happens crossing the non-interacting point c = 0. When entering the attractive phase, bound states are expected to be produced and we proposed an educated exact ansatz for this process. However, aside the low-density limit, its validity is extremely challenging to be tested due to the absence of independent analytical or numerical computations.

In this work, we address this protocol in the classical NLS equation: the motivations are multifaceted. First of all, by comparing with large-scale classical simulations we can benchmark the analytical predictions, exploring corrections to adiabaticity and due to integrability breaking corrections, always present in experimental setups. Secondly, the in se-per-se interest in the thermodynamic description of integrable models is fueled by the advances of GHD studies. In contrast with quantum models, where excitations are of fermionic or bosonic nature, classical systems feature either solitonic or radiative modes, obeying Maxwell–Boltzmann or Rayleigh–Jeans statistics. Probably, the simplest example of classical integrability is the hard rods gas [63], namely an ensemble of hard-core particles of finite length obeying a (dressed) Maxwell–Boltzmann statistics, thus having solitonic nature. Solitonic modes also appear in the Toda chain [64, 65] and Landau–Lifschitz spin model [66–68]. In contrast, the thermodynamics of the relativistic sinh-Gordon model [69, 70] and the repulsive NLS [71] is built upon radiative modes. The GHD of an integrable discretization of the NLS has been recently addressed [72] and, unexpectedly, the excitation spectrum has been shown to be described by solitonic modes in analogy with the Toda chain, but in sharp contrast with the naive expectation dragged from the continuous case. The soliton-radiation correspondence in the continuum limit remains yet to be clarified.

Within the attractive phase, the continuous NLS is known [27] to support both radiative and solitonic modes, hence its thermodynamics can be naively expected to be described by a mixture of the two. Surprisingly, we show that this is not the case: by accessing the classical GHD through a semiclassical limit of the quantum model, we recover only solitonic modes. Benchmarks with ab initio classical simulations confirm the correctness of the procedure.

It is worth mentioning that far before the advent of GHD, the soliton picture (see reference [73] for a recent review) provides a hydrodynamic framework to study transport in classical integrable models bearing solitonic excitations. Indeed, the GHD equations describing inhomogeneous states, but with homogeneous and time-independent dynamics, reduce to the soliton gas picture when applicable. Nonetheless, passing through the semiclassical limit of the quantum model bears many advantages, such as describing the effect of interaction changes. Most importantly for our purposes, with this strategy we recover for free the correct excitations that describe the thermodynamics and hydrodynamics of the classical model. Before going through the details of our results, we wish to clarify a possible source of confusion: several studies have been dedicated to the so-called semiclassical (or weakly dispersive) limit of the NLS and its hydrodynamic form (see e.g. reference [74]). In our approach, we consider the semiclassical limit of the quantum Bose gas to access a description of the NLS valid in arbitrary regimes, albeit in the GHD limit.

The paper is organized as follows. Section 2 reviews the solution of the quantum Bose gas, its hydrodynamics and slow interaction changes from the repulsive to the attractive phase. Section 3 presents the classical NLS, discussing the difficulties of obtaining a hydrodynamic description from purely classical methods and thus building the GHD from semiclassical limits of the quantum model. Section 4 benchmarks the hydrodynamic treatment against ab initio Monte Carlo simulations. Finally in section 5 we gather our conclusions followed by two short appendixes describing the numerical methods.

In this section, we summarize the key points of the exact solution of the Bose gas and its hydrodynamic treatment, which we use as a stepping stone to access the classical NLS. We start by reviewing the microscopic solution by means of the coordinate Bethe ansatz and discussing the strings representing bound-states in the attractive phase. We then turn to the thermodynamic limit and conclude with a short summary of GHD and its application to slow drives from a repulsive to an attractive interaction strength.

2.1. The many-body wavefunction

Already in 1963 Lieb and Liniger provided the exact solution of the homogeneous V(x) = 0 1D Bose gas [32], which is also referred to as LL model. Integrable models are characterized by the presence of an extensive number of local conservation laws which greatly affects the dynamics. Indeed, the Hilbert space can be understood as a collection of particles undergoing pairwise elastic scatterings. The rapidities λ_i represent the asymptotic momenta of the particles after a scattering event. The asymptotic states are common eigenvectors $\vert {\left\{{\lambda }_{i}\right\}}_{i=1}^{N}\rangle$ of the Hamiltonian and all the charges, which act additively on the quasiparticles due to their locality

$$\begin{equation}\hat{\mathcal{Q}}\vert {\left\{{\lambda }_{i}\right\}}_{i=1}^{N}\rangle =\left(\sum\limits _{i=1}^{N}\enspace q({\lambda }_{i})\right)\vert {\left\{{\lambda }_{i}\right\}}_{i=1}^{N}\rangle ,\end{equation} \tag{ 3 }$$

where q(λ) is called charge eigenvalue. Notable examples are the number of particles q(λ) → 1, the energy q(λ) → (λ) = λ² and the momentum q(λ) → p(λ) = λ.

The many-body wavefunction is explicitly computed via the coordinate Bethe ansatz [32] and reads

$$\begin{equation}{\Psi}({\left\{{x}_{i}\right\}}_{i=1}^{N};{\left\{{\lambda }_{i}\right\}}_{i=1}^{N})=\sum\limits _{\mathcal{P}}\enspace A(\mathcal{P})\prod\limits _{j=1}^{N}{\text{e}}^{\text{i}{\lambda }_{{\mathcal{P}}_{j}}{x}_{j}},\quad {x}_{1}\leqslant {x}_{2}\leqslant \dots \leqslant {x}_{N}.\end{equation} \tag{ 4 }$$

Above, $\mathcal{P}$ are the permutations of rapidities and the coefficients $A(\mathcal{P})$ encode the factorized scattering of the particles. In particular, exchanging two rapidities is equivalent to a scattering process, hence the state picks up a phase. By defining Π_j,j+1 as the permutation that exchanges the rapidities in position j and j + 1, it holds that $A({{\Pi}}_{j,j+1}\mathcal{P})=S({\lambda }_{{\mathcal{P}}_{j}}-{\lambda }_{{\mathcal{P}}_{j+1}})A(\mathcal{P})$, with

$$\begin{equation}S(\lambda )\equiv {\text{e}}^{\text{i}{\Theta}(\lambda )}\equiv \frac{\lambda +\mathrm{i}c}{\lambda -\mathrm{i}c},\end{equation} \tag{ 5 }$$

where S(λ) and Θ(λ) are the two-body scattering matrix and phase, respectively. Notice, that for any non zero interaction, one has S(0) = −1: as a consequence, the wavefunction (4) vanishes if two or more identical rapidities are considered, enforcing a Pauli exclusion principle on the set ${\left\{{\lambda }_{i}\right\}}_{i=1}^{N}$. As we will see, the scattering phase fully encodes the interacting nature of the problem and it is a key ingredient in its thermodynamic description. The first step in this direction is to consider a finite system of length L: periodic boundary conditions quantize the allowed rapidities through the Bethe equations

$$\begin{equation}{\mathrm{e}}^{\mathrm{i}L{\lambda }_{j}}=\prod\limits _{k\ne j}S({\lambda }_{j}-{\lambda }_{k}),\quad j,k=1,\dots ,N.\end{equation} \tag{ 6 }$$

Notice that in the non-interacting limit S(λ) = 1 the standard quantization conditions of non-interacting models are recovered ${\mathrm{e}}^{\mathrm{i}L{\lambda }_{j}}=1$. In the presence of interactions, the highly non-linear nature of the Bethe equations make their solution extremely challenging.

In the zero density limit (L → ∞ and N fixed) the nature of the solutions can be classified and strongly depends on the sign of the interaction strength. Within the repulsive phase c > 0, only real rapidities λ_j are admitted, which describe asymptotically free particles. On the contrary, in the attractive regime c < 0 bound states are possible and are described by complex rapidities. In fact, complex λ_j with the proper imaginary part are associated with a wavefunction that exponentially vanishes for large particle separation (4).

At the level of the Bethe equation (6), a complex rapidity causes ${\text{e}}^{\text{i}{\lambda }_{j}L}$ to either diverge or vanish when L → +∞. The zero or divergence must be counterbalanced by a similar behavior in the product of the scattering matrices and this locks the rapidity differences to either the zeroes or poles of S(λ), giving rise to regular patterns in the complex plane which are called strings. Rapidities belonging to the same string have the same real part, but are regularly shifted along the imaginary direction ${\left\{\lambda -i\frac{c}{2}(j+1-2a)\right\}}_{a=1}^{j}$ where j ∈ {1, 2, ...} is the length of the string. Physically, a string represents a bound state of j particles traveling together with momentum λ. The string hypothesis [34] assumes that in the true thermodynamic limit (L → ∞ with N/L kept fixed), the strings can still be used to classify the state. Hence, by grouping together the constituents of each string, one can write the Bethe–Takahashi equations

$$\begin{equation}{\text{e}}^{\text{i}L{p}_{j}({\lambda }_{j}^{\alpha })}=\prod\limits _{({k}^{\prime },{\alpha }^{\prime })\ne (j,\alpha )}\enspace {S}_{j{j}^{\prime }}({\lambda }_{j}^{\alpha }-{\lambda }_{{j}^{\prime }}^{{\alpha }^{\prime }}),\end{equation} \tag{ 7 }$$

which constrain the real center of mass ${\lambda }_{j}^{\alpha }$ of the strings where j specifies the string species. The string scattering matrix ${S}_{j{j}^{\prime }}(\lambda )={\text{e}}^{\text{i}{{\Theta}}_{j{j}^{\prime }}(\lambda )}$ is simply defined by the product of the scattering matrices of all the constituents of the two scattered strings. In terms of the scattering phase, it reads

$$\begin{equation}{{\Theta}}_{jk}(\lambda )=(1-{\delta }_{jk}){\theta }_{\vert j-k\vert }(\lambda )+2{\theta }_{\vert j-k\vert +2}(\lambda )+\cdots +2{\theta }_{j+k-2}(\lambda )+{\theta }_{j+k}(\lambda )\end{equation} \tag{ 8 }$$

$$\begin{equation}\text{and}\enspace {\theta }_{j}(\lambda )=-2\enspace \mathrm{arctan}\frac{2\lambda }{cj}.\end{equation} \tag{ 9 }$$

Similarly, one defines the charge eigenvalues of a given string by summing the contributions of its constituents ${q}_{j}(\lambda )={\sum }_{a=1}^{j}\enspace q\left(\right.\lambda -i\frac{c}{2}\left(\right.j+1-2a\left.\right)\left.\right)$. In particular, in the case of energy and momentum one has

$$\begin{equation}{{\epsilon}}_{j}(\lambda )=j{\lambda }^{2}-\frac{{c}^{2}}{12}j({j}^{2}-1),\qquad {p}_{j}(\lambda )=j\lambda .\end{equation} \tag{ 10 }$$

Particularly relevant is the energy eigenvalue from which one can readily extract important information on the stability of the attractive phase. In fact, the attractive ground state is a giant bound state of all its constituents with energy E_g = −c² N(N² − 1)/12. Remaining in the zero momentum sector, the first excited state is obtained by removing one particle from the molecule and placing it infinitely far apart, hence E_e = −c²(N − 1)[(N − 1)² − 1]/12. Thus, the energy difference E_e − E_g diverges in the thermodynamic limit, signaling the impossibility of exciting the system for any finite temperature.

2.2. The Bose gas in the thermodynamic limit

The thermodynamic limit is attained by sending the system size to infinity L → ∞ while retaining a fixed density n = N/L. In this perspective, solving the Bethe equations by analytical or numerical methods appears a daunting task which, however, can be avoided in the thermodynamic limit. In the spirit of the thermodynamic Bethe ansatz (TBA) [34], the exact details of a microstate defined by a set of discrete rapidities $\vert {\left\{{\lambda }_{i}\right\}}_{i=1}^{N}\rangle$ are no longer of importance. Instead, we pass over to a coarse-grained description of states in terms of the so-called root densities ρ_j (λ). For the sake of a more compact discussion, we consider the repulsive and attractive case at the same time, hence we introduce a set of root densities {ρ_j (λ)}, one for each string species: the repulsive phase is obtained by retaining only strings with j = 1. Concretely, Ldλρ_j (λ) counts how many solutions to the Bethe–Takahashi equations of the jth string specie lay within the interval [λ, λ + dλ], thus equation (3) in the thermodynamic limit becomes

$$\begin{equation}\langle {\left\{{\lambda }_{i}\right\}}_{i=1}^{N}\vert \hat{\mathcal{Q}}\vert {\left\{{\lambda }_{i}\right\}}_{i=1}^{N}\rangle \to L\sum\limits _{j}\int \mathrm{d}\lambda \enspace {q}_{j}(\lambda ){\rho }_{j}(\lambda ).\end{equation} \tag{ 11 }$$

The Bethe–Takahashi equation (7) in logarithmic form constrain the allowed set of rapidities $2\pi {\mathcal{I}}_{j,\alpha }=L{p}_{j}({\lambda }_{j}^{\alpha })-{\sum }_{({k}^{\prime },{\alpha }^{\prime })\ne (j,\alpha )}{{\Theta}}_{j,{j}^{\prime }}({\lambda }_{j}^{\alpha }-{\lambda }_{{j}^{\prime }}^{{\alpha }^{\prime }})$, with ${\mathcal{I}}_{j,\alpha }$ an integer. As a consequence, one introduces the total root density ${\rho }_{j}^{t}(\lambda )\geqslant {\rho }_{j}(\lambda )$ as the maximum allowed occupancy and ${\vartheta }_{j}(\lambda )={\rho }_{j}(\lambda )/{\rho }_{j}^{t}(\lambda )\leqslant 1$ is the filling fraction. Due to the non-linear nature of the Bethe–Takahashi equations, ${\rho }_{j}^{t}(\lambda )$ is explicitly state-dependent and is defined as the Jacobian of the change of coordinates between the rapidities and the set of integers ${\mathcal{I}}_{j,\alpha }$. In the thermodynamic limit, this is compactly expressed as $2\pi {\rho }_{j}^{t}(\lambda )={({\partial }_{\lambda }{p}_{j})}^{\mathrm{d}r}$, where the dressing operation is introduced for an arbitrary test function ${\tau }_{j}(\lambda )\to {\tau }_{j}^{\mathrm{d}r}(\lambda )$ as

$$\begin{equation}{\tau }_{j}^{\mathrm{d}r}(\lambda )={\tau }_{j}(\lambda )-\sum\limits _{k}\int \frac{\mathrm{d}{\lambda }^{\prime }}{2\pi }{\varphi }_{jk}(\lambda -{\lambda }^{\prime }){\vartheta }_{k}({\lambda }^{\prime }){\tau }_{k}^{\mathrm{d}r}({\lambda }^{\prime }).\end{equation} \tag{ 12 }$$

Here, the kernel φ_jk (λ) is the derivative of the scattering phase φ_jk (λ) = ∂_λ Θ_jk (λ). The dressing operation is ubiquitous in TBA calculations, as we will see. Due to the coarse-grained nature of the root density representation, many states share the same macroscopic description: the Yang–Yang entropy

$$\begin{equation}\mathcal{S}[{\left\{{\rho }_{j}\right\}}_{j=1}^{\infty }]=\int \frac{\mathrm{d}\lambda }{2\pi }\sum\limits _{j}\;{\rho }_{j}^{t}(\lambda )\left(-{\vartheta }_{j}\enspace \mathrm{ln}\enspace {\vartheta }_{j}-(1-{\vartheta }_{j})\mathrm{ln}(1-{\vartheta }_{j})\right)\end{equation} \tag{ 13 }$$

counts the degeneracy of the macrostate identified by a given root density, in the sense that the number of microscopic states sharing the same root density scales as $\sim \mathrm{exp}(L\mathcal{S}[{\left\{{\rho }_{j}\right\}}_{j=1}^{\infty }])$. Notice that $\mathcal{S}$ is nothing else than the entropy of 'fermionic particles' with occupancy ϑ_j (λ) and total phase space distribution ${\rho }_{j}^{t}(\lambda )$: the fermionic nature of the Yang–Yang entropy is due to the Pauli principle exclusions enforced on the rapidities. The entropic contribution is crucial when aiming to describe an ensemble of states, both in equilibrium and non-equilibrium setups.

Equilibrium thermodynamics —one of the earliest successes of integrability has been its ability to exactly describe systems at thermal equilibrium, despite the presence of strong interactions. Let us focus on a thermal state ${\text{e}}^{-\beta (\hat{H}-\mu \hat{N})}$ with β and μ the inverse temperature and chemical potential, respectively. The grand canonical partition function $\mathcal{Z}$ is obtained by summing over all microscopic states. In the thermodynamic limit, one can coarse grain the microscopic summation by replacing it with a path integral over the root densities

$$\begin{equation}\mathcal{Z}=\sum\limits _{N,{\left\{{\lambda }_{i}\right\}}_{i=1}^{N}}{\mathrm{e}}^{-\sum\limits _{i=1}^{N}\left(\right.\beta ({\epsilon}({\lambda }_{i})-\mu )}\to \int \mathcal{D}{\rho }_{j}\enspace {\mathrm{e}}^{L\left[\mathcal{S}[{\rho }_{j}]-\sum\limits _{j}\int \mathrm{d}\lambda \beta ({{\epsilon}}_{j}(\lambda )-\mu ){\rho }_{j}(\lambda )\right]}.\end{equation} \tag{ 14 }$$

In the thermodynamic limit L → ∞, the path integral localizes on the saddle point solution obtained by minimizing the free energy $F=-\mathcal{S}[\left\{{\rho }_{j}\right\}]+{\sum }_{j}\int \mathrm{d}\lambda \enspace \beta [{{\epsilon}}_{j}(\lambda )-\mu ]{\rho }_{j}(\lambda )$. This minimization procedure determines a unique choice of root densities that fully describes the thermal state. Within the attractive regime, the ground state's instability is also present at finite temperatures, therefore the resulting equations are stable only within the repulsive regime and read [34]

$$\begin{equation}\varepsilon (\lambda )=\beta [{\epsilon}(\lambda )-\mu ]-\int \frac{\mathrm{d}{\lambda }^{\prime }}{2\pi }\varphi (\lambda -{\lambda }^{\prime })\mathrm{log}(1+{\mathrm{e}}^{-\varepsilon (\lambda )}),\end{equation} \tag{ 15 }$$

where one conveniently parameterizes the filling fractions through the effective energy ɛ

$$\begin{equation}\vartheta (\lambda )=\frac{1}{1+{\mathrm{e}}^{\varepsilon (\lambda )}}.\end{equation} \tag{ 16 }$$

Generalized Gibbs ensembles —besides playing a pivotal role at equilibrium, the root densities fully characterize several non-equilibrium protocols. Let us focus on the paradigmatic example of the homogeneous quantum quench [75]: here, the system is initialized in a certain state (or ensemble) that is not an eigenstate of the Hamiltonian. Then, the state evolves under unitary evolution and eventually (local) relaxation is attained. The presence of infinitely many local conservation laws greatly constrains the dynamics of integrable models and hinders thermalization. Indeed, the steady state reached by integrable models after a quantum quench is usually not described by thermal ensembles, but rather by GGE [50] ${\mathrm{e}}^{-{\sum }_{i}\;{\beta }_{i}{\hat{\mathcal{Q}}}_{i}}$, where in addition to the Hamiltonian, all the relevant charges must be taken into account. Within the thermodynamic limit, GGEs can be described with minor changes of the thermal analysis, by modifying the source term in the TBA integral equation (15): however, there are two major obstructions to this program. The first concerns the charges that need to be included in the GGE: in continuum models such as the Bose gas, the expectation values of ultra local charges often feature ultraviolet divergences [76, 77], which must be properly regularized [78–80]. In addition, also in the absence of ultraviolet divergences as for example in lattice systems, listing all the charges needed in the GGE is challenging and quasi-local conservation laws must be considered [81, 82]. Fortunately, the explicit construction of the charges is not needed for the thermodynamic description of the GGE. Indeed, GGEs are understood to be in one-to-one correspondence with the root density description [83]: given a set of root densities, a GGE of the form ${\mathrm{e}}^{-{\sum }_{i}\;{\beta }_{i}{\hat{\mathcal{Q}}}_{i}}$ exists from which they can be derived.

The second challenge to be addressed is determining the final GGE from the initial conditions: a major breakthrough in this direction has been the quench action approach [84, 85] which extracts the root densities from the knowledge of the overlaps between the initial state and the post quench eigenstates. Thanks to this approach, several quench protocols in integrable models have been exactly solved, but its applicability is limited to those initial states where the aforementioned overlap can be computed [86].

In the case of the Bose gas, the list of these special states hosts a unique element, namely the non-interacting ground state (i.e. the homogeneous Bose–Einstein condensate). The sudden quench from the non-interacting ground state to the repulsive case is among the first quenches in interacting integrable models ever solved [87] and it has been lately generalized to the case of attractive interactions [60, 61], showing bound state formation. However, despite the appeal of such an exact solution, the instability of the one-dimensional BEC against thermal fluctuations makes the realization of this protocol in real-life experiments difficult. In contrast to sudden quantum quenches, slow time-dependent protocols allows for a much greater flexibility in the choice of the initial state: these setups fall within the realm of hydrodynamics, which we now discuss.

2.3. Generalized hydrodynamics and bound state formation

When a many-body system is subjected to slow and smoothly varying external perturbations, it can be rightfully expected to locally evolve through quasi-stationary states. This is the onset of hydrodynamics: by invoking a separation of time scales, one assumes the system is locally and instantaneously described by a homogeneous stationary state, which then evolves on macroscopic times through proper continuity equations. For those models that are evolving through a locally integrable Hamiltonian, stationary states are described by GGEs and the hydrodynamics of GGEs is now known as GHD [36, 37]. Since its introduction in 2016, it has provided tremendous advances in understanding 1D systems out of equilibrium. Here, we focus on the hydrodynamic equations at the Euler scale, i.e. we retain only the first order in a gradient expansion: diffusive corrections are known [88, 89] and play a pivotal role in the late-time thermalization due to the trap-induced integrability breaking corrections [90]. The thermalization time scale is proportional to the square of the sloshing period and can be neglected for very smooth traps. At the Eulerian scale, the hydrodynamics is governed by the following continuity equation

$$\begin{equation}{\partial }_{t}{\rho }_{j}+{\partial }_{x}({v}_{j}^{\text{eff}}{\rho }_{j})+{\partial }_{\lambda }({F}_{j}^{\text{eff}}{\rho }_{j})=0,\end{equation} \tag{ 17 }$$

where the space-time dependent root density describes the local GGE and we suppress the explicit dependence on rapidities, space and time for the sake of a more compact notation. We discuss the attractive and repulsive cases again at once, the latter being obtained by retaining only the first string species ρ_j=1 in the forthcoming expressions. The continuity equation describes the phase space density of particles moving with velocity ${v}_{j}^{\text{eff}}$ and experiencing a force ${F}_{j}^{\text{eff}}$: due to the strongly interacting nature of the system. The velocity and force are suitably renormalized by the presence of the surrounding particles. The form ${v}_{j}^{\text{eff}}={({\partial }_{\lambda }{{\epsilon}}_{j})}^{\mathrm{d}r}/{({\partial }_{\lambda }{p}_{j})}^{\mathrm{d}r}$, where the dressing operation is defined in equation (12), has been proposed in the original papers and then later rigorously proven in references [91, 92] (however, see reference [93] for an exception). The effective force undergoes an analogous dressing operation ${F}_{j}^{\text{eff}}={f}_{j}^{\mathrm{d}r}/{({\partial }_{\lambda }{p}_{j})}^{\mathrm{d}r}$, where f_j is the bare force, which captures the effects of traps [45] and interaction changes [46]

$$\begin{equation}{f}_{j}(\lambda )=-j{\partial }_{x}V+{\partial }_{t}c\sum\limits _{{j}^{\prime }}\int \mathrm{d}{\lambda }^{\prime }\enspace {\partial }_{c}{{\Theta}}_{j{j}^{\prime }}(\lambda -{\lambda }^{\prime }){\rho }_{{j}^{\prime }}({\lambda }^{\prime }).\end{equation} \tag{ 18 }$$

Above, we assume for simplicity that the interaction changes in time, but not in space. The hydrodynamic equation allows for evolving the root density within the repulsive or the attractive phase, but does not provide a recipe to connect the two. For example, let us initialize the system in a given state (e.g. thermal) within the repulsive phase c > 0 and then gently tune the interaction to smaller values until c = 0⁺ is met. Here, the evolution with the repulsive GHD stops and continues from c = 0⁻ within the attractive regime with the initial conditions being determined by the endpoint of the repulsive evolution. This boundary condition is not contained within the GHD equations presented so far, which must be supplemented by further considerations: in reference [62], the three of us proposed an analytical ansatz to accomplish this task (see also [42, 94]). Our reasoning is based on entropic arguments and on the fact that at c = 0 bound states have purely real rapidities and are indistinguishable from unbound particles with the same rapidity. This last statement is expressed by the continuity equation

$$\begin{equation}{\left.\rho (\lambda )\right\vert }_{c={0}^{+}}=\sum\limits _{j}\enspace j{\left.{\rho }_{j}(\lambda )\right\vert }_{c={0}^{-}},\end{equation} \tag{ 19 }$$

which, besides being diagonal in coordinate space (whose dependence we omit) is also diagonal in rapidity space. In the case of an interaction change from attractive to repulsive, the right hand-side is known and uniquely determines the repulsive root density: physically, a bound state of j particles traveling with velocity ${v}_{j}^{\text{eff}}(\lambda )$ melts in j unbound particles traveling with the same velocity ${v}^{\text{eff}}({\lambda }^{\prime })={v}_{j}^{\text{eff}}(\lambda )$. At the non-interacting point, it is easy to check that the last identity implies λ = λ'. On the other hand, when the interaction is tuned from repulsive to attractive, only particles with the same rapidity can form a bound state, since particles with different velocity would drift apart before having the chance to bind. Notice that in the more interesting case when the attractive phase is approached from the repulsive one, equation (19) does not uniquely fix {ρ_j (λ)} and further considerations must be made. If bound states at c = 0⁻ are indistinguishable from unbound particles, any choice of the bound states population is equally probable: therefore, the thermodynamic limit will select the most probable configuration, obtained by entropy maximization. Maximizing the Yang–Yang entropy (13) at the free point under the constraint of particle conservation (19) uniquely determines the attractive root densities through

$$\begin{equation}{\varepsilon }_{j}(\lambda )=j\omega (\lambda )+\sum\limits _{k}(2\enspace \mathrm{min}(j,k)-{\delta }_{jk})\mathrm{log}(1+{\mathrm{e}}^{-{\varepsilon }_{k}(\lambda )}),\end{equation} \tag{ 20 }$$

where ω(λ) is the rapidity-dependent Lagrange multiplier whose value must be tuned to enforce equation (19). The effective energy ɛ_j parameterizes the state through ${\vartheta }_{j}(\lambda )={(1+{\mathrm{e}}^{{\varepsilon }_{j}(\lambda )})}^{-1}$.

Similar equations already appeared in the literature [95, 96] in the context of the non-interacting limit of the sine-Gordon (SG) field theory and can be analytically solved ⁴

$$\begin{equation}{\vartheta }_{j}(\lambda )=1-\frac{\mathrm{sinh}(j\omega (\lambda )/2)\mathrm{sinh}((j+2)\omega (\lambda )/2)}{{\mathrm{sinh}}^{2}((j+1)\omega (\lambda )/2)}.\end{equation} \tag{ 21 }$$

As a last step, the Lagrange multiplier ω(λ) can be connected with the repulsive root density through $\rho (\lambda )=\frac{1}{2\pi }\frac{1}{{\text{e}}^{\omega (\lambda )}-1}$, obtaining an explicit solution for the filling fraction.

The GHD equations together with the ansatz (20) allow one to access the attractive phase in a controlled manner from a large variety of initial states. Once c is gently tuned to negative values, the gas enters a highly excited phase featuring bound states of finite size, whose stability is protected by integrability. In contrast to equilibrium thermal states, GGEs obtained by solving (20) are stable and the bound state population is exponentially decaying in their size ϑ_j (λ) ∼ e^−jω(λ). Of course, the value of ω(λ) > 0 is fixed imposing equation (19): it decreases monotonously as $\rho (\lambda ){\vert }_{c={0}^{+}}$, approaching zero when $\rho (\lambda ){\vert }_{c={0}^{+}}$ diverges and leading to instabilities. Our ansatz is rooted on physical arguments and well-motivated thermodynamic considerations, but ab initio checks are extremely difficult. Analytical consistency checks are available only in the low density limit, while numerical benchmarks are a daunting task. The state-of-the-art tensor network simulations struggle to capture the dynamics of highly excited and entangled systems, especially in the case where the model is continuous and of bosonic nature (which results in an infinite local Hilbert space). In fact, the appeal of GHD is to provide trustful predictions in a regime beyond any other method. As we will see in the next section, the same hydrodynamic treatment together with the bound states production ansatz can be transferred to the classical world: there, efficient numerical simulations will confirm the correctness of our ansatz.

After having reviewed the quantum Bose gas and its hydrodynamic description, we now revert to the classical world and introduce the NLS equation equation (2). The equation of motion can be derived from a classical Hamiltonian

$$\begin{equation}H=\int \mathrm{d}x\enspace \vert {\partial }_{x}\phi {\vert }^{2}+\varkappa \vert \phi {\vert }^{4}+V(x)\vert \phi {\vert }^{2},\end{equation} \tag{ 22 }$$

by enforcing the Poisson brackets {ϕ(x), ϕ*(y)} = iδ(x − y). From the comparison between the Hamiltonian of the 1D quantum model (1) with that of the classical NLS, the similarities between the two appear self-evident. Indeed, the homogeneous NLS (V(x) = 0) is a well-known integrable partial differential equation [26, 27]: it features infinitely many local conservation laws affecting the dynamics. In the absence of a Hilbert space, there is not a classical analogue of the coordinate Bethe ansatz and the NLS is studied with the inverse scattering method (ISM). In contrast, the quantum model can also be studied by means of the quantum ISM [97], but its discussion is beyond the scope of our work. The key point is noticing that the NLS equation emerges as a compatibility condition of an auxiliary linear problem

$$\begin{equation}{\partial }_{x}F(t,x)={U}_{\lambda }(t,x)F(t,x)\enspace ,\qquad {\partial }_{t}F(t,x)={V}_{\lambda }(t,x)F(t,x),\end{equation} \tag{ 23 }$$

where F is an auxiliary two-dimensional vector and U_λ , V_λ are 2 × 2 matrices with an explicit dependence on ϕ

$$\begin{equation}{U}_{\lambda }=\sqrt{\varkappa }\left({\phi }^{\ast }{\sigma }^{+}+\phi {\sigma }^{-}\right)+\frac{\lambda }{2i}{\sigma }^{z}\end{equation} \tag{ 24a }$$

$$\begin{equation}{V}_{\lambda }=\left(i\varkappa \vert \phi {\vert }^{2}+{\lambda }^{2}\frac{i}{2}\right){\sigma }^{z}-i\sqrt{\varkappa }\left({\partial }_{x}\phi {\sigma }^{+}-{\partial }_{x}\phi {\sigma }^{-}\right)-\lambda \sqrt{\varkappa }({\phi }^{\ast }{\sigma }^{+}+\phi {\sigma }^{-})\end{equation} \tag{ 24b }$$

with σ^x,y,z the standard Pauli matrices and σ^± = (σ^x ± iσ^y )/2.

The consistency of the two differential equation (23) requires the zero curvature condition ∂_t U_λ − ∂_x V_λ + [U_λ , V_λ ] = 0 to hold, the latter is λ-independent and it is satisfied if and only if the Schrӧdinger equation (2) holds. Eventually, we are interested in the emergent thermodynamic properties in setups with finite particle number and energy densities, i.e. where L⁻¹∫dx|ϕ|² remains constant in the thermodynamic limit. However, it is convenient to start our discussion within the zero density limit. An extensive discussion can be found in references [26, 27].

The whole line monodromy matrix —in this case, the system is assumed to be infinite with the condition that the wavefunction ϕ(t, x) vanishes sufficiently fast (e.g. exponentially) for |x| → +∞. Let us define the transfer matrix as the spatial propagator of equation (23)

$$\begin{equation}T(t,x,y)=\text{Pexp}{\int }_{x}^{y}\mathrm{d}z\enspace {U}_{\lambda }(t,z),\end{equation} \tag{ 25 }$$

where Pexp is the path-ordered exponential. Thanks to the zero curvature condition, the transfer matrix obeys a simple equation of motion

$$\begin{equation}{\partial }_{t}{T}_{\lambda }(t,x,y)={V}_{\lambda }(t,x){T}_{\lambda }(t,x,y)-{T}_{\lambda }(t,x,y){V}_{\lambda }(t,y).\end{equation} \tag{ 26 }$$

Particularly important is the role of the monodromy matrix, i.e. the extension of the transfer matrix to the whole line

$$\begin{equation}{\mathcal{T}}_{\lambda }(t)=\underset{x\to \infty }{\mathrm{lim}}\;{\mathrm{e}}^{\frac{\lambda x}{2i}{\sigma }^{z}}{T}_{\lambda }(t,-x,x){\mathrm{e}}^{-\frac{\lambda x}{2i}{\sigma }^{z}}.\end{equation} \tag{ 27 }$$

Above, the additional oscillating terms are introduced in view of the asymptotic behavior ${\mathrm{lim}}_{\vert x\vert \to \infty }{U}_{\lambda }(t,x)=\lambda {\sigma }^{z}/(2i)$. Due to the vanishing boundary conditions, one simply has ${\mathrm{lim}}_{\vert x\vert \to +\infty }{V}_{\lambda }(t,x)=i{\sigma }^{z}{\lambda }^{2}/2$. From equation (26) it immediately follows ${\partial }_{t}\enspace \mathrm{Tr}[{\mathcal{T}}_{\lambda }]=0$ for any spectral parameter. As a consequence, an infinite λ-dependent family of conservation laws is found. Further insight is gained from the matrix elements of the monodromy matrix, which can be parameterized as

$$\begin{equation}{\mathcal{T}}_{\lambda }=\left(\begin{matrix}\hfill a(\lambda )\hfill & \hfill \mathrm{sign}(\varkappa ){b}^{\ast }(\lambda )\hfill \\ \hfill b(\lambda )\hfill & \hfill {a}^{\ast }(\lambda )\hfill \end{matrix}\right),\end{equation} \tag{ 28 }$$

where a(λ), b(λ) are entire functions in the λ-complex plane and play the role of action-angle variables of the integrable system. They in fact obey a very simple equation of motion

$$\begin{equation}{\partial }_{t}a(\lambda )=0\enspace ,\qquad {\partial }_{t}b(\lambda )=-i{\lambda }^{2}b(\lambda ).\end{equation} \tag{ 29 }$$

Hence, a(λ) is itself a conserved quantity. In principle, the non-linear dynamics of the wavefunction ϕ(t, x) is solved, albeit in an implicit way. In practice, one should compute the monodromy matrix from the initial profile, then evolve the scattering parameters according to equation (29) and finally use them to infer back the scattering potential and thus ϕ. The last step is what is referred to as 'inverse scattering' and passes through the solution of the Gelfand–Levitan–Marchenko linear integral equation. In this respect, the analytical properties of a(λ) are crucial. Indeed, a(λ) can be analytically continued in the complex and a(λ) → 1 when λ is large, leading to the representation

$$\begin{equation}a(\lambda )=\prod\limits _{a=1}\frac{\lambda -{\lambda }_{a}}{\lambda +{\lambda }_{a}^{\ast }}\enspace \mathrm{exp}\left[{\int }_{-\infty }^{\infty }\frac{\mathrm{d}{\lambda }^{\prime }}{\pi i}\frac{\mathrm{log}\vert a({\lambda }^{\prime })\vert }{{\lambda }^{\prime }-\lambda -i{0}^{+}}\right].\end{equation} \tag{ 30 }$$

Above, $\mathfrak{I}({\lambda }_{a}) > 0$ and the complex zeroes are associated with solitons, i.e. non-dispersive solutions of the NLS traveling at constant velocity. Solitons are present only for attractive interactions ϰ and are instead absent in the repulsive case. In contrast with solitons, radiative modes (i.e. dispersive solutions) are present for both the repulsive and attractive case. Equation (30) has a suggestive interpretation: let us consider the repulsive phase where solitons are absent and expand log(λ) for large values of λ

$$\begin{equation}\enspace \mathrm{log}\enspace a(\lambda )=i\sum\limits _{n=0}^{\infty }\enspace {\lambda }^{-n-1}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}{\lambda }^{\prime }}{\pi }{({\lambda }^{\prime })}^{n}\enspace \mathrm{log}\vert a({\lambda }^{\prime })\vert .\end{equation} \tag{ 31 }$$

On the one hand, a(λ) is an integral of motion for any λ in view of equation (29), hence each term ${\mathcal{Q}}_{n}\equiv {\int }_{-\infty }^{\infty }\frac{\mathrm{d}{\lambda }^{\prime }}{\varkappa }{({\lambda }^{\prime })}^{n}\enspace \mathrm{log}\vert a({\lambda }^{\prime })\vert$ is separately conserved. In addition, it is possible to show that ${\mathcal{Q}}_{n}$ are local conserved charges [26, 69]. On the other hand, a comparison with the expectation value of the quantum charges in the thermodynamic limit (11) suggests that log |a(λ)| plays the role of a classical root density [69, 98]. However, a proper analysis requires to consider the system in a finite volume, which we now briefly discuss.

The periodic monodromy matrix —when periodic boundary conditions are enforced, the monodromy matrix is defined as the transfer matrix across the finite system ${\mathcal{T}}_{\lambda }^{\text{PBC}}(t)={T}_{\lambda }(t,-L/2,L/2)$. Equation (26) still holds and by virtue of PBC V_λ (t, −L/2) = V_λ (t, L/2), hence the trace of the monodromy matrix is still a conserved quantity. There are two main differences to the previous case: first, not all choices of the spectral parameter are acceptable. In fact, by imposing periodic boundary conditions only values of λ such that ${\mathcal{T}}_{\lambda }^{\text{PBC}}(t)=\pm \text{Id}$ are allowed. Secondly, in contrast with the fast decaying case, only the trace of the monodromy matrix is conserved, while the diagonal entries are not. However, this is due to boundary terms that are unimportant in the thermodynamic limit if one focuses only on bulk properties: following the argument of references [69, 98], let us imagine a configuration of the field ϕ initially having support on the restricted interval [−L/2 + ℓ, L/2 − ℓ]. In this case, the monodromy matrix on the whole line coincides with the PBC case and has the same analytic properties of its coefficient. In addition, a(λ) is conserved until the evolving field does not reach the boundaries. This suggests that in the thermodynamic limit one can identify ${[{\mathcal{T}}_{\lambda }^{\text{PBC}}]}_{1,1}\simeq {\text{e}}^{-\text{i}\lambda L/2}a(\lambda )$ with a(λ) having the same analytic properties of the whole line problem. Focusing on the repulsive case, this observation led to the identification of the classical analogue of the Bethe equation (6): indeed, by imposing PBC $\mathrm{log}{[{\mathcal{T}}_{{\lambda }_{j}}^{\text{PBC}}]}_{1,1}=\mathrm{log}(\pm 1)=\pi j$ and using the representation equation (30) (where solitons are absent since ϰ > 0), one obtains the set of quantization conditions

$$\begin{equation}{\lambda }_{j}-{\unicode{10767}}_{-\infty }^{\infty }\mathrm{d}{\lambda }^{\prime }\frac{2\varkappa }{{\lambda }_{j}-{\lambda }^{\prime }}\frac{\mathrm{log}\vert a({\lambda }^{\prime })\vert }{\varkappa \pi L}=2\pi j{L}^{-1},\end{equation} \tag{ 32 }$$

with ⨏ being the principal value regularization of the integral. These equations have strong analogies with the quantum Bethe equations in the logarithmic form (6), where the spectral parameter is identified with the rapidities and the summation over rapidities has already been replaced with an integral over the root density. In particular, the classical root density ρ(λ) and scattering phase Θ(λ) can be identified as [69]

$$\begin{equation}\rho (\lambda ){\leftrightarrow}\frac{1}{\pi \varkappa L}\enspace \mathrm{log}\vert a(\lambda )\vert \enspace ,\qquad {\Theta}(\lambda ){\leftrightarrow}\frac{2\varkappa }{\lambda }.\end{equation} \tag{ 33 }$$

The proportionality factor between ρ(λ) and log |a(λ)| is determined comparing the local charges obtained from equation (31) with the standard normalization of particle density and Hamiltonian. This correspondence leads to a full characterization of the post-quench GGE in the repulsive phase [71] in terms of the initial data. In principle, one could attempt the same strategy in the attractive case as well, but this path appears unpractical or, at least, very challenging. First, solitonic modes complicate the analytical properties of the scattering coefficient and must be properly taken into account. Strictly speaking, a sharp distinction between solitons and radiation exists only in the whole line problem and is blurred out in the periodic case: nevertheless, we still refer to the complex zeroes of a(λ) as solitons. Secondarily, in writing equation (32) one implicitly assumes only real values of the spectral parameter λ_j . However, in the attractive case one expects complex solutions as well, in analogy with the quantum case and their classification further complicates the approach. Fortunately, one can avoid a solution of the classical problem and indirectly access its thermodynamic description through a proper semiclassical limit of the well-controlled quantum case.

The NLS as the semiclassical limit of the quantum Bose gas —despite the intrinsic quantumness of nature, our daily life experience is based on classical physics: in the large mode occupation limit, quantum systems are trustfully described by classical equations of motion. Semiclassical methods have a long lasting success story which is not possible to properly credit here, therefore we restrict our attention to the interplay between semiclassical limits and integrability, emphasizing the connection between the NLS and the quantum Bose gas. At equilibrium, the quantum–classical correspondence has been studied in early days [95, 99], but it was brought back into focus through recent interest in non-equilibrium problems. In reference [69] the idea of studying classical GGEs by taking the semiclassical limit of quantum theories has been put forward, which has been later successfully applied to sudden quench protocols in the NLS [71]. Following the advent of GHD, the hydrodynamics of classical integrable models has been successfully derived from their quantum relatives [67, 70], in particular the hydrodynamics of the repulsive NLS has been proven to be extremely fruitful in benchmarking GHD predictions against ab initio numerical simulations. Here, we will follow a similar route and tackle the NLS in its attractive phase ϰ < 0, which becomes accessible in view of our previous results in the attractive quantum Bose gas [62]. The semiclassical limit is most conveniently accessed by introducing a small parameter h which plays the role of Planck's constant. The large occupation limit requires the number of particles to diverge N ∝ h⁻¹, but a finite limit is attained only if the quantum interaction is rescaled to small values c = hϰ. This is readily seen from a path integral representation of the quantum propagator

$$\begin{align}\hfill \langle {\psi }_{j}\vert {\mathrm{e}}^{-iT\hat{H}}\vert {\psi }_{i}\rangle & ={\left.\int \mathcal{D}\psi \enspace {\text{e}}^{\text{i}{\int }_{0}^{T}\mathrm{d}{t}^{\prime }\int \mathrm{d}x\enspace \frac{i}{2}({\partial }_{t}{\psi }^{{\dagger}}\psi -{\psi }^{{\dagger}}{\partial }_{t}\psi )-\vert {\partial }_{x}\psi {\vert }^{2}-c\vert \psi {\vert }^{4}}\right\vert }_{\psi (x,T)={\psi }_{f}(x),\psi (x,0)={\psi }_{i}(x)}=\hfill \\ \hfill & \hspace{-9pt}\stackrel{\psi =\phi /\sqrt{h}}{=}{\left.\int \mathcal{D}\phi \enspace {\mathrm{e}}^{{h}^{-1}i{\int }_{0}^{T}\mathrm{d}{t}^{\prime }\int \mathrm{d}x\enspace \frac{i}{2}({\partial }_{t}{\phi }^{{\dagger}}\phi -{\phi }^{{\dagger}}{\partial }_{t}\phi )-\vert {\partial }_{x}\phi {\vert }^{2}-\varkappa \vert \phi {\vert }^{4}}\right\vert }_{\phi (x,T)={h}^{-1/2}{\psi }_{f}(x),\phi (x,0)={h}^{-1/2}{\psi }_{i}(x)}.\hfill \end{align} \tag{ 34 }$$

In the h → 0 limit, the path integral localizes on the saddle point, enforcing the classical equations of motion. The whole analysis is then carried out on thermal ensembles and GGEs. In particular, the limit of large occupation modes must be imposed for each conserved charge, leading to the scaling

$$\begin{equation}\langle {\mathcal{Q}}^{\text{quantum}}\rangle ={h}^{-1}\langle {\mathcal{Q}}^{\text{classical}}\rangle .\end{equation} \tag{ 35 }$$

In the next section, we briefly review the semiclassical limit of the quantum TBA and GHD in the repulsive case before moving on with the attractive interactions.

3.1. The semi-classical limit of the repulsive LL

In this and in the forthcoming sections, it is convenient to introduce an extra label to tell apart quantum and classical quantities. Therefore, we add the subscript 'q' to variables labelling quantum quantities. For the repulsive phase, we follow references [69, 71] (see also reference [100] for an alternative derivation). As we anticipated, the classical model features a single root density describing radiative modes. Let us rewrite the large occupancy scaling of the charges (35) using the root density representation

$$\begin{equation}\int \mathrm{d}\lambda \enspace {q}_{\mathrm{q}}(\lambda ){\rho }_{\mathrm{q}}(\lambda )={h}^{-1}\int \mathrm{d}\lambda \enspace q(\lambda )\rho (\lambda ).\end{equation} \tag{ 36 }$$

The charge eigenvalues of the local charges such as the density or the energy do not bear any explicit interaction dependence. Hence, one can identify the quantum eigenvalues with the classical ones q_q(λ) = q(λ): assuming this identification holds for a complete set of charges, imposing equation (36) and at the leading order one is led to the scaling ρ_q(λ) = h⁻¹ ρ(λ), with possible further corrections in h. Indeed, this quantum–classical correspondence turns out to give a finite limit of the whole TBA with the already anticipated classical scattering phase (33). The classical scattering phase is recovered by the quantum one as

$$\begin{equation}{\Theta}(\lambda )=\underset{h\to 0}{\mathrm{lim}}\frac{1}{h}\left({{\Theta}}_{\mathrm{q}}(\lambda )+\pi \enspace \mathrm{sign}(\lambda )\right).\end{equation} \tag{ 37 }$$

In order to keep the standard definition of the total root density, the classical root density is consequentially obtained as ${\rho }^{t}(\lambda )={\mathrm{lim}}_{h\to 0}({\rho }_{\mathrm{q}}^{t}(\lambda )-{\rho }_{\mathrm{q}}(\lambda ))$. With this caveat, the classical dressing is defined in the usual manner employing the classical scattering phase and the identification 2πρ^t (λ) = (1)^dr (λ) still holds: consistently, the limit of the whole GHD results in the same equations as equation (17) where the classical scattering phase replaces the quantum one. An important difference arises in the integral equations describing thermal states and this is deeply connected to the different statistical properties of the classical NLS when compared with the quantum Bose gas, encoded in equation (37). Indeed, the Pauli principle exclusion on the quantum Bethe roots naturally leads to fermionic statistics (13), while the radiative modes of the repulsive NLS are expected to follow a Rayleigh–Jeans distribution, which can be thought of as the large occupation limit of bosonic particles. In fact, equation (37) removes the fermionic component of the scattering phase before sending h → 0: this is strictly connected with the possibility of reformulating the whole quantum TBA in terms of bosonic degrees of freedom [101]. By taking the semi-classical limit of the Yang–Yang entropy (13), one finds

$$\begin{equation}{\mathcal{S}}_{\mathrm{q}}\simeq \int \frac{\mathrm{d}\lambda }{2\pi }\left\{{\rho }^{t}(\lambda )\mathrm{log}\enspace \vartheta (\lambda )+(1-\mathrm{log}\enspace h){\rho }^{t}(\lambda )\right\},\end{equation} \tag{ 38 }$$

that still has a divergent term ∝log h. However, by explicitly using the definition of ρ^t it is immediate to see that ∫dλρ^t (λ) is a state-independent, since

$$\begin{equation}\int \mathrm{d}\lambda \enspace \left\{{\rho }^{t}(\lambda )-\frac{1}{2\pi }\right\}=\int \frac{\mathrm{d}\lambda }{2\pi }{\partial }_{\lambda }{\Theta}(\lambda -{\lambda }^{\prime })\rho ({\lambda }^{\prime })=0.\end{equation} \tag{ 39 }$$

Hence, it can be neglected, leading to the definition of the classical Yang–Yang entropy within the repulsive phase

$$\begin{equation}\mathcal{S}\equiv \int \frac{\mathrm{d}\lambda }{2\pi }{\rho }^{t}(\lambda )\mathrm{log}\enspace \vartheta (\lambda ).\end{equation} \tag{ 40 }$$

Just as in the quantum case, thermal states are obtained by minimizing the free energy, leading to a Rayleigh–Jeans distribution of the mode occupation ϑ(λ) = 1/ɛ(λ)

$$\begin{equation}\varepsilon (\lambda )=\beta \left[{\epsilon}(\lambda )-\mu \right]-{\unicode{10767}}_{-\infty }^{\infty }\frac{\mathrm{d}{\lambda }^{\prime }}{2\pi }{\partial }_{\lambda }{\Theta}(\lambda -{\lambda }^{\prime })\mathrm{log}\enspace \varepsilon ({\lambda }^{\prime }).\end{equation} \tag{ 41 }$$

Thanks to the equipartition theorem, non-interacting waves on thermal states are distributed according to the inverse of the energy: equation (41) generalizes this expectation to the presence of interactions. Equation (41) matches the natural expectation inferred by the inverse scattering solution, from which only radiative modes of a single species are expected. This must be contrasted with the attractive phase we now discuss.

3.2. The attractive phase

As we have already discussed, tackling the thermodynamics and hydrodynamics of the attractive phase with the inverse scattering solution appears challenging. Nevertheless, the semiclassical limit offers a direct path to our goal, whose validity has to be ultimately benchmark against numerical simulations. The first point to be addressed is the correspondence between the quantum and the classical root densities: we start by noticing that in the semiclassical limit the energy quantization of the strings is expected to merge in a continuum. In order to understand the correct scaling, let us focus on the expectation value of a given quantum charge

$$\begin{align}\hfill {L}^{-1}\langle \hat{\mathcal{Q}}\rangle & =\sum\limits _{j=1}^{\infty }\int \mathrm{d}\lambda \enspace {q}_{j}(\lambda ){\rho }_{\mathrm{q};j}(\lambda )\hfill \\ \hfill & =\sum\limits _{j=1}^{\infty }\int \mathrm{d}\lambda \enspace \sum\limits _{a=1}^{j}\enspace q(\lambda -ih(j+1-2a)\varkappa /2){\rho }_{\mathrm{q};j}(\lambda )\sim {h}^{-1},\hfill \end{align} \tag{ 42 }$$

where we explicated the charge eigenvalue using the string hypothesis as well as the correspondence between the quantum and classical interaction strength c = hϰ. By requiring the same scaling $\sim {h}^{-1}$ for all the charges and assuming their completeness, the same scaling must be independently enforced on each rapidity λ. Furthermore, the imaginary shifts in the charge eigenvalue are equispaced on a vanishing step $\sim h$. This suggests to replace the summation over the quantum strings with an integral over a continuum variable s = hj. Therefore we get

$$\begin{align}\hfill & \sum\limits _{j=1}^{\infty }\enspace \sum\limits _{a=1}^{j}\enspace q(\lambda -ih(j+1-2a)\varkappa /2){\rho }_{\mathrm{q};j}(\lambda )\hfill \\ \hfill & \qquad \to {\int }_{0}^{\infty }\mathrm{d}s{\int }_{0}^{s}\mathrm{d}{s}^{\prime }q(\lambda -i(s-2{s}^{\prime })\varkappa /2){h}^{-2}{\rho }_{\mathrm{q};s/h}.\hfill \end{align} \tag{ 43 }$$

Imposing the scaling $\sim {h}^{-1}$ of the whole expression suggests the following quantum–classical correspondence

$$\begin{equation}{\rho }_{\mathrm{q};j}(\lambda )=h{\rho }_{hj}(\lambda ).\end{equation} \tag{ 44 }$$

Following this argument, the attractive phase accommodates for the semiclassical limit in a completely different manner when compared with the repulsive case: instead of macroscopically populating a single root density, it fills a divergent number of particle species (i.e. the bound states), but each of them with a population of order $\mathcal{O}(h)$. This is a priori not entirely obvious from equation (43), since the same scaling could have been obtained with a single non-vanishing and divergent summand. However, equation (43) is not the only requirement to be enforced and equation (44) guarantees a finite limit of all the TBA equations, as we now discuss.

Before turning to the total root density, we can first deduce the semiclassical limit of the quantum scattering phase (8), which explicitly reads

$$\begin{align}\hfill {{\Theta}}_{\mathrm{q};jk}(\lambda )& =({\delta }_{jk}-1)2\enspace \mathrm{arctan}\left(\frac{2\lambda }{\varkappa h\vert j-k\vert }\right)-2\enspace \mathrm{arctan}\left(\frac{2\lambda }{\varkappa h(j+k)}\right)\hfill \\ \hfill & \quad -\sum\limits _{a=1}^{\mathrm{max}(j,k)-1}\enspace 4\enspace \mathrm{arctan}\left(\frac{2\lambda }{\varkappa h(\vert j-k\vert +2a)}\right).\hfill \end{align} \tag{ 45 }$$

Notice, that by taking the limit of large strings necessary for attaining the semiclassical limit, the summation diverges as $\sim {h}^{-1}$, while the first two terms remain $\mathcal{O}({h}^{0})$. Therefore, we define the classical scattering phase as ${{\Theta}}_{s{s}^{\prime }}(\lambda )={\mathrm{lim}}_{h\to 0}\;h{{\Theta}}_{\mathrm{q};s/h\enspace \enspace {s}^{\prime }/h}(\lambda )$: the first two terms in equation (45) do not contribute in this limit

$$\begin{equation}{{\Theta}}_{s{s}^{\prime }}(\lambda )=-2{\int }_{\vert s-{s}^{\prime }\vert }^{s+{s}^{\prime }}\mathrm{d}\sigma \enspace \mathrm{arctan}\left(\frac{2\lambda }{\varkappa \sigma }\right).\end{equation} \tag{ 46 }$$

Using the classical scattering phase, one can write the total root density, the dressing equations and all the GHD equations. Of course, they coincide with the semiclassical limit of the quantum hydrodynamics. In particular, the quantum total root density scales as ${\rho }_{\mathrm{q};j}^{t}(\lambda )=h{\rho }_{hj}^{t}(\lambda )$, with the classical total root density defined as the dressed derivative of the momentum $2\pi {\rho }_{s}^{t}(\lambda )={({\partial }_{\lambda }{p}_{s}(\lambda ))}^{\mathrm{d}r}$. The classical dressing for an arbitrary test function ${\tau }_{s}\to {\tau }_{s}^{\mathrm{d}r}$ is defined through the usual integral equation

$$\begin{equation}{\tau }_{s}^{\mathrm{d}r}(\lambda )={\tau }_{s}(\lambda )-{\int }_{0}^{\infty }\mathrm{d}{s}^{\prime }\int \mathrm{d}{\lambda }^{\prime }{\varphi }_{s{s}^{\prime }}(\lambda -{\lambda }^{\prime }){\vartheta }_{{s}^{\prime }}({\lambda }^{\prime }){\tau }_{{s}^{\prime }}^{\mathrm{d}r}({\lambda }^{\prime }),\end{equation} \tag{ 47 }$$

with ${\vartheta }_{s}(\lambda )={\rho }_{s}(\lambda )/{\rho }_{s}^{t}(\lambda )$ the classical filling function and φ_ss'(λ) = ∂_λ Θ_ss'(λ). It is worth to mention that φ_ss'(λ) = ∂_λ Θ_ss'(λ) can be explicitly written as

$$\begin{equation}{\varphi }_{s{s}^{\prime }}(\lambda )=\frac{2}{\varkappa }\enspace \mathrm{ln}\frac{{(s+{s}^{\prime })}^{2}{\varkappa }^{2}+4{\lambda }^{2}}{{(s-{s}^{\prime })}^{2}{\varkappa }^{2}+4{\lambda }^{2}}.\end{equation} \tag{ 48 }$$

The scattering kernel φ_ss' (48) is nothing else than the post-scattering displacement of two colliding classical solitons [102] and has recently appeared in connection with superdiffusion in the XXZ spin chain [39, 67], where the exotic transport has been attributed to excitations of semiclassical nature (see also references [68, 103, 104]). In the semiclassical limit, the XXZ spin chain is described by the Landau–Lifschitz model, which is gauge-equivalent to the attractive NLS [27, 105], hence the identical scattering kernels. The energy and momentum eigenvalues are immediately obtained as ${p}_{s}(\lambda )={\mathrm{lim}}_{h\to 0}\;h{p}_{\mathrm{q};s/h}(\lambda )$ and ${{\epsilon}}_{s}(\lambda )={\mathrm{lim}}_{h\to 0}\;h{{\epsilon}}_{\mathrm{q};s/h}(\lambda )$ and read

$$\begin{equation}{p}_{s}(\lambda )=s\lambda \enspace ,\qquad {{\epsilon}}_{s}(\lambda )=s{\lambda }^{2}-\frac{{\varkappa }^{2}}{12}{s}^{3}.\end{equation} \tag{ 49 }$$

Finally, the GHD equations have the usual form ${\partial }_{t}{\rho }_{s}+{\partial }_{x}({v}_{s}^{\text{eff}}{\rho }_{s})+{\partial }_{\lambda }({F}_{s}^{\text{eff}}{\rho }_{s})=0$, with the following effective quantities ${v}_{s}^{\text{eff}}={({\partial }_{\lambda }{{\epsilon}}_{s})}^{\mathrm{d}r}/{({\partial }_{\lambda }{p}_{s})}^{\mathrm{d}r}$, ${F}_{s}^{\text{eff}}={f}_{s}^{\mathrm{d}r}/{({\partial }_{\lambda }{p}_{s})}^{\mathrm{d}r}$ and the bare force

$$\begin{equation}{f}_{s}(\lambda )=-s{\partial }_{x}V+{\partial }_{t}\varkappa {\int }_{0}^{\infty }\mathrm{d}{s}^{\prime }\enspace {\partial }_{\varkappa }{{\Theta}}_{s,{s}^{\prime }}(\lambda -{\lambda }^{\prime }){\rho }_{{s}^{\prime }}({\lambda }^{\prime }).\end{equation} \tag{ 50 }$$

We now turn to the semi-classical limit of the Yang–Yang entropy, leading to the definition of the classical entropy within the attractive phase. On the practical side, this will allow us to access the attractive regime by means of slow interaction changes. Moreover, it will unveil the nature of the excitations described by ρ_s (λ).

By using the quantum–classical correspondence in equation (13), one finds

$$\begin{equation}\mathcal{S}=\underset{\delta (h) > 0{}}{\overset{+\infty }{\int }}\mathrm{d}s\int \mathrm{d}\lambda \enspace {\rho }_{s}^{t}(\lambda )\left\{{\vartheta }_{s}(\lambda )[1-\mathrm{log}\enspace {\vartheta }_{s}(\lambda )]-2\enspace \mathrm{log}\enspace (h)\enspace {\vartheta }_{s}(\lambda )\right\},\end{equation} \tag{ 51 }$$

where a $\sim \mathrm{log}\enspace h$ term survives and the role of δ(h) will be clarified soon. In contrast with the repulsive phase, where the entropy was pointing out the radiative nature of the excitations, equation (51) is associated with solitons. Indeed, consider N indistinguishable classical particles, hence the associated phase space is 1/N! and the entropy reads $\mathcal{S}=\mathrm{log}(1/N!)\sim N(1-\mathrm{log}\enspace N)$, in striking similarity with equation (51). We stress the absence of explicit radiative modes, as one could have naively expected from the ISM. Let us now analyze the role of the $\sim \mathrm{log}\enspace h$ term: a similar correction appeared in the repulsive case as well, but, as we commented, its contribution was state independent and could be omitted. This is not the case in equation (51) and, in order for the h → 0 limit to give a finite equation of state through entropy maximization, the $\sim \mathrm{log}\enspace h$ term should be counterbalanced by another singularity. This is the role covered by δ(h): indeed, as we now discuss, ϑ_s develops a power law singularity for small s and δ(h) must be properly tuned to obtain a finite expression. From the semiclassical limit perspective, we are summing over string indexes j = {1, 2, 3, ...} hence s ∈ {h, 2h, 3h, ...}: therefore, a cutoff δ(h) ∝ h is expected. We notice that a similar log h term appeared also in the semiclassical TBA of the XXZ chain [67], further strengthening the connection with the NLS. In that case, the role of the Planck constant was played by an external magnetic field. In the following, we determine δ(h) through a bootstrap approach, asking the maximum entropy approach to guide us towards a finite (and explicitly h-independent) formulation of the bound state production.

Crossing ϰ = 0 and soliton production —when crossing over from the repulsive to the attractive phase, the classical problem can be tackled in the same spirit as the quantum Bose gas. One can proceed in two different but equivalent ways, either by taking directly the semiclassical limit of the exact quantum solution (21) or through maximization of the classical entropy. We start by addressing the first case, then we turn back to the classical entropy maximization and discuss the role of the log(h) term appearing in equation (51). The semiclassical limit of equation (21) is extremely simple to be taken: by enforcing the h-scaling discussed in the previous section, one reaches the simple expression

$$\begin{equation}{\vartheta }_{s}(\lambda )=\frac{1}{{(2\pi \rho (\lambda ))}^{2}}\frac{1}{4\enspace {\mathrm{sinh}}^{2}\left(\frac{s}{4\pi \rho (\lambda )}\right)}.\end{equation} \tag{ 52 }$$

Notice that the filling fraction diverges for small s as ϑ_s (λ) ∼ 1/s², regardless of the value of ρ(λ). This divergence is the reason behind the appearance of log(h) in the classical entropy, as we now discuss. By enforcing that the two phases should be indistinguishable when approaching the non-interacting point, one reaches the semiclassical analogue of (19)

$$\begin{equation}{\left.\rho (\lambda )\right\vert }_{\varkappa ={0}^{+}}={\int }_{0}^{\infty }\mathrm{d}s\enspace s{\left.{\rho }_{s}(\lambda )\right\vert }_{\varkappa ={0}^{-}}.\end{equation} \tag{ 53 }$$

This continuity equation does not suffice to determine the population of each string, which is in turn determined by the constrained maximization of the classical entropy. In the limit ϰ → 0⁻, the kernel of the dressing equation becomes diagonal in rapidity space

$$\begin{equation}\underset{\varkappa \to {0}^{-}}{\mathrm{lim}}{\varphi }_{s,{s}^{\prime }}(\lambda )=\underset{\varkappa \to {0}^{-}}{\mathrm{lim}}\frac{2}{\varkappa }\enspace \mathrm{ln}\frac{{\varkappa }^{2}{\left(s+{s}^{\prime }\right)}^{2}+4{\lambda }^{2}}{{\varkappa }^{2}{\left(s-{s}^{\prime }\right)}^{2}+4{\lambda }^{2}}=-4\pi \mathrm{min}(s,{s}^{\prime })\delta (\lambda ).\end{equation} \tag{ 54 }$$

As a consequence, the definition of the total root density is also diagonal in rapidity space and the entropy can be maximized for each rapidity independently leading to the equation

$$\begin{equation}{\varepsilon }_{s}(\lambda )=s\omega (\lambda )+2\enspace \mathrm{log}\enspace h+{\int }_{\delta (h)}^{\infty }\mathrm{d}{s}^{\prime }2\enspace \mathrm{min}(s,{s}^{\prime }){\mathrm{e}}^{-{\varepsilon }_{{s}^{\prime }}(\lambda )},\end{equation} \tag{ 55 }$$

where we used the natural parameterization ${\vartheta }_{s}(\lambda )={\mathrm{e}}^{-{\varepsilon }_{s}(\lambda )}$. The λ-dependent Lagrange multiplier ω(λ) is of course determined by enforcing equation (53). These integral equation can develop singular solutions for s → 0: this is precisely what happens in equation (52) and the divergence must be singled out before taking the h → 0 limit. In view of equation (52), it is useful to introduce the reparameterization

$$\begin{equation}{\varepsilon }_{s}=2\enspace \mathrm{log}\enspace s+s{\bar{\varepsilon }}_{s},\end{equation} \tag{ 56 }$$

where ${\bar{\varepsilon }}_{s}$ remains finite as s → 0. Then, the domain [δ, +∞] in the integral of equation (55) is split into two parts [δ, C] ∪ [C, +∞]. The constant C will be later sent to zero and it is meant to be small enough such that $s{\bar{\varepsilon }}_{s}\simeq 0$ for s ∈ [δ, C]. In equation (55) we now focus on s ∈ [δ, C] and use equation (56)

$$\begin{align}\hfill 2\enspace \mathrm{log}\enspace s+s{\bar{\varepsilon }}_{s}& =s\omega (\lambda )+2\enspace \mathrm{log}\enspace h+{\int }_{\delta (h)}^{C}\mathrm{d}{s}^{\prime }2\frac{\mathrm{min}(s,{s}^{\prime })}{{s}^{\prime 2}}+{\int }_{\delta (h)}^{C}\mathrm{d}{s}^{\prime }2\enspace \mathrm{min}(s,{s}^{\prime })\frac{{\mathrm{e}}^{-{s}^{\prime }{\bar{\varepsilon }}_{{s}^{\prime }}(\lambda )}-1}{{s}^{\prime 2}}\hfill \\ \hfill & \quad +s{\int }_{C}^{\infty }\mathrm{d}{s}^{\prime }2\frac{{\mathrm{e}}^{-{s}^{\prime }{\bar{\varepsilon }}_{{s}^{\prime }}(\lambda )}}{{s}^{\prime 2}}.\hfill \end{align} \tag{ 57 }$$

In the above, the last integral is obviously finite. The first integral on the second row is not only finite, but it vanishes as C → δ → 0. Hence, the singularity is entirely due to the integral in the first row ${\int }_{\delta (h)}^{C}\mathrm{d}{s}^{\prime }2\frac{\mathrm{min}(s,{s}^{\prime })}{{s}^{\prime 2}}=2\enspace \mathrm{log}(s/\delta )+2s\left(\frac{1}{s}-\frac{1}{C}\right)=2\enspace \mathrm{log}(s/\delta (h))+2+\mathcal{O}(s)$. By matching the singular and $\mathcal{O}({s}^{0})$ terms in equation (57), we are lead to log δ(h) = log h + 1. The next natural question concerns the effects of the singular filling ϑ_s on the behavior of the total root density and, in general, of dressed quantities. A similar analysis of the s ∼ 0 behavior can be performed, showing that the total root density is not only finite, but it also quadratically vanishes (see next subsection) as s → 0: this is crucial, since it guarantees the root density ${\rho }_{s}(\lambda )={\rho }_{s}^{t}{\vartheta }_{s}(\lambda )$ to remain finite, thus giving finite expectation values of the charges. Using the identity log δ = log h + 1, one can obtain a finite integral equation for $\bar{\varepsilon }$ and explicitly take the h → 0 limit. Indeed, by using ${\int }_{\delta (h)}^{\infty }\mathrm{d}{s}^{\prime }2\frac{\mathrm{min}(s,{s}^{\prime })}{{s}^{\prime 2}}=2\enspace \mathrm{log}(s/\delta (h))+2$, we can rewrite equation (55) as

$$\begin{equation}{\bar{\varepsilon }}_{s}(\lambda )=\omega (\lambda )+{\int }_{\delta (h)}^{\infty }\mathrm{d}{s}^{\prime }\frac{2\enspace \mathrm{min}(s,{s}^{\prime })}{s}\frac{{\mathrm{e}}^{-{s}^{\prime }{\bar{\varepsilon }}_{{s}^{\prime }}(\lambda )}-1}{{s}^{\prime 2}}.\end{equation} \tag{ 58 }$$

The integrand is now finite in zero and the limit δ(h) → 0 can be safely taken. As in the quantum case, this integral equation admits an exact analytical solution which, with the identification $\rho (\lambda )=\frac{1}{2\pi \omega (\lambda )}$, is equivalent to the result (52) derived as the semiclassical limit of the quantum filling.

In figure 1 we compare the semiclassical scaling with the quantum result for small values of h and study the convergence. In the case of the filling fractions, we compare the analytical solutions while the correspondent root densities are obtained by numerically solving the dressing equations. Overall, we experience that a large number of quantum strings are needed to reproduce the classical result, especially in the case of the filling function. In contrast, the classical dressing equations can be efficiently discretized with a far lower number of points (in this example, excellent convergence is already obtained with $\sim 20$ points). This stresses the need of solving directly the classical GHD rather than numerically considering the semiclassical limit of the quantum case.

Zoom In Zoom Out Reset image size

A non-singular reparameterization —the analysis of the singularities we provided suggests to redefine the total root density and filling by singling out the small s behavior. The same procedure is also extremely useful for computing quantities appearing in the GHD equations. Hence, we introduce the following definitions that are also valid for finite interactions ϰ < 0.

$$\begin{equation}{\bar{\vartheta }}_{s}(\lambda )\equiv {s}^{2}{\vartheta }_{s}(\lambda )\enspace ,\qquad {\bar{\rho }}_{s}^{t}(\lambda )\equiv {s}^{-2}{\rho }_{s}^{t}(\lambda ).\end{equation} \tag{ 59 }$$

Furthermore, we define the regularized dressing on a test function ${t}_{s}(\lambda )\to {t}_{s}^{\mathbf{d}\mathbf{r}}(\lambda )$ as the solution of the following integral equation

$$\begin{equation}s{t}_{s}^{\mathbf{d}\mathbf{r}}(\lambda )={t}_{s}(\lambda )-{\int }_{0}^{\infty }\mathrm{d}{s}^{\prime }\int \mathrm{d}{\lambda }^{\prime }{s}^{-1}{\varphi }_{s{s}^{\prime }}(\lambda -{\lambda }^{\prime }){\bar{\vartheta }}_{{s}^{\prime }}({\lambda }^{\prime }){t}_{{s}^{\prime }}^{\mathbf{d}\mathbf{r}}({\lambda }^{\prime }).\end{equation} \tag{ 60 }$$

Notice that the limit ${\mathrm{lim}}_{s\to 0}\;{s}^{-1}\enspace {\varphi }_{s{s}^{\prime }}(\lambda )$ is finite and non-vanishing for any interactions. As a consequence, if the source term does not have singularities ${t}_{s}(\lambda )\sim \mathcal{O}({s}^{0})$ as s → 0, then ${t}_{s}^{\mathbf{d}\mathbf{r}}(\lambda )$ is expected to be finite as well, which is easily numerically checked. The definition equation (60) leads to the following correspondence between the regularized and standard dressing (47)

$$\begin{equation}{\tau }_{s}^{\mathrm{d}r}(\lambda )={s}^{2}{({s}^{-1}{\tau }_{s})}^{\mathbf{d}\mathbf{r}}\end{equation} \tag{ 61 }$$

In particular, this implies ${\bar{\rho }}_{s}^{t}(\lambda )={(1/(2\pi ))}_{s}^{\mathbf{d}\mathbf{r}}(\lambda )$, hence ${\bar{\rho }}_{s}^{t}(\lambda )$ remains finite as s → 0, while ${\rho }_{s}^{t}(\lambda )$ quadratically vanishes, as anticipated. Due to the $\sim s$ behavior of the energy and bare force, ${({s}^{-1}{\partial }_{\lambda }{{\epsilon}}_{s})}^{\mathbf{d}\mathbf{r}}$ and ${({s}^{-1}{f}_{s})}^{\mathbf{d}\mathbf{r}}$ are finite as well. Hence, the effective velocity and force can now be conveniently computed in terms of a ratio of non-singular functions.

One of the major advantages that classical field theories have in contrast to quantum ones is the possibility of numerically accessing arbitrarily long times, hence we are in the perfect position to test our hydrodynamic predictions and quantify the corrections. In particular, classical simulations can be performed through Monte Carlo methods: the initial field configurations are randomly sampled from a known probability distribution. Then, the initial conditions are deterministically evolved through the equations of motion and the desired observables are computed. Finally, the average over the initial conditions is taken. Thermal states within the repulsive phase are probably the most natural choice for the initial conditions and can be efficiently simulated by means of metropolis-hasting methods [106–108]. However, as we will discuss in the following, other choices of initial GGEs are also plausible and easy to simulate. We perform the standard microscopic numerical simulations following references [71, 90]. For the interested reader, a short summary is reported in appendix A. When numerically solving the GHD equations, it turns out that the root density is not the most convenient variable. Already in the seminal papers [36, 37] the GHD equations have been rewritten in terms of the filling fraction where convective form is recovered. Within the attractive phase, the GHD equations in terms of the filling fraction are

$$\begin{equation}{\partial }_{t}{\bar{\vartheta }}_{s}+{v}_{s}^{\text{eff}}{\partial }_{x}{\bar{\vartheta }}_{s}+{F}_{s}^{\text{eff}}{\partial }_{\lambda }{\bar{\vartheta }}_{s}=0,\end{equation} \tag{ 62 }$$

and an analogous expression is reached within the repulsive phase, where of course the string index is absent. Notice that, importantly, it is convenient to use the rescaled filling fraction defined in equation (59) to ensure the regularity of the expression for s → 0. Similarly, in the attractive phase the effective velocity and force are computed using the modified definition of the dressing introduced in equations (60) and (61). Showing the equivalence between the GHD equation (62) and the canonical version in terms of the root density requires to manipulate the integral equations defining the dressing operation. We choose to not report here these lengthy, but standard calculations: the interested reader can find them in the seminal papers [36, 37] (see also reference [46] for the presence of forces) and straightforwardly apply them to our classical case as well. Once the GHD equations are written in the form equation (62), they can be seen as infinitesimal translations in the (x, λ) plane and solved by the method of characteristics [46, 109]: we provide a short summary in appendix B. In principle, the GHD equation (62) can be used to study interaction changes in fully inhomogeneous settings, that are for example given by trapping potentials. In practice, we experience that the need of discretizing the continuum set of strings of the attractive phase places a major hurdle on the numerical simulations. For this reason, and since we are mostly interested in observing the bound state production due to time-dependent interaction changes, we focus here on the homogeneous case. Nevertheless, a more careful optimization of the GHD code would make the fully inhomogeneous case accessible as well. Notice that, within the homogeneous case, the ∂_t ϰ term appearing in the effective force can be singled out and absorbed with a change of variables. Therefore, the GHD evolution can be parameterized in terms of the interaction only. In the following, we use this convention. The comparison between hydrodynamics and microscopic simulations passes through the expectation value of observables. As we discussed, in integrable models the expectation values of conserved quantities are the most easy ones to be accessed, leading to the density of particles and energy as the natural choices. However, the particle density is trivially conserved throughout the protocol and is thus not informative at all. On the other hand, the expectation value of the energy density diverges on classical thermal states: this is due to the well known ultraviolet catastrophe of black body radiation, which led to the discovery of quantum mechanics. Indeed, at very high momenta one can discard the role of interactions in equation (41) obtaining ρ(k) ∼ (β(k))⁻¹ and causing the energy density to diverge ∫dk(k)ρ(k) → ∞. For this reason, we focus on the expectation value of the correlated density ⟨|ϕ(x)|⁴⟩, which is finite and undergoes a non-trivial evolution. In the quantum case, its expectation value can be computed on arbitrary GGEs by means of a generalization of the Hellmann–Feynmann theorem (see e.g. [62]). From this point onwards, the semiclassical limit is readily taken leading to the classical expressions. In particular, one has

$$\begin{equation}\langle \vert \phi {\vert }^{4}\rangle =\int \frac{\mathrm{d}\lambda }{2\pi }2\lambda \vartheta (\lambda ){f}^{\mathrm{d}\mathrm{r}}(\lambda ),\hspace{24.0pt}\varkappa > 0\end{equation} \tag{ 63 }$$

$$\begin{equation}\langle \vert \phi {\vert }^{4}\rangle =\int \mathrm{d}s\mathrm{d}\lambda \left\{-{\rho }_{s}(\lambda )\frac{\varkappa }{6}{s}^{3}+\frac{2s\lambda }{2\pi }{\vartheta }_{s}(\lambda ){f}_{s}^{\mathrm{d}\mathrm{r}}(\lambda )\right\}\hspace{24.0pt}\varkappa < 0.\end{equation} \tag{ 64 }$$

In the second line we can use the non-singular parameterization we introduced by noticing ${\vartheta }_{s}(\lambda ){f}_{s}^{\mathrm{d}\mathrm{r}}(\lambda )={\bar{\vartheta }}_{s}(\lambda ){({s}^{-1}{f}_{s})}^{\mathbf{d}\mathbf{r}}$. The amount of information of the evolving state captured by the correlated density is of course limited and a better characterization of the state though further measurement is desirable. Within the repulsive case, close analytical expression for the exact probability distribution of the density |ϕ|² have been recently derived [71]. Consequentially, all the momenta ⟨|ϕ|²ⁿ ⟩ can be computed on arbitrary GGEs (see also references [110, 111] for the quantum case). Moreover, the root density itself could be determined from the microscopic simulations thanks to the correspondence (33). However, to the best of our knowledge, analogous results are not available within the attractive realm, hence we turn our attention to different observables in the form of two-point correlation functions. Indeed, GHD allows for an exact computation of two-point correlation functions $\langle \mathcal{O}(t,x)\mathcal{O}({t}^{\prime },y)\rangle$ on the Eulerian scale (i.e. leading order in the large space and time separation limit), provided the expectation value of $\langle \mathcal{O}\rangle$ on arbitrary GGEs is known. This result has first been obtained for the time independent and homogeneous GGEs [112] and was tested for the first time on the classical Sinh Gordon model [70]. Later, expressions for inhomogeneous and time dependent GGEs became available [113] (see also [38]). For simplicity, here we focus on the two-point correlation of the density on constant GGEs

$$\begin{equation}t\langle \mathcal{O}(t,x)\mathcal{O}{(0,0)\rangle }_{\text{c}}=\int \mathrm{d}s\enspace \mathrm{d}\lambda \enspace \enspace \delta \left(x/t-{v}_{s}^{\text{eff}}(\lambda )\right){\rho }_{s}(\lambda ){\sigma }_{s}(\lambda ){[{V}_{s}^{\mathcal{O}}(\lambda )]}^{2},\end{equation} \tag{ 65 }$$

where ⟨...⟩_c stands for the connected part of the correlator and the integration over the string variable is obviously omitted in the repulsive case ϰ > 0. Equation (65) holds for large space-time separations and after averaging over mesoscopic fluid cells. Notice the scaling form of the right hand term as a function of the ray ζ ≡ x/t. The function ${V}^{\mathcal{O}}$ carries the information of the observable of interest and it is defined as ${\rho }_{s}^{t}(\lambda ){V}_{s}^{\mathcal{O}}(\lambda )=\frac{\delta \langle \mathcal{O}\rangle }{\delta {\vartheta }_{s}(\lambda )}$. Notably, in the case where $\mathcal{O}$ is the local density of a conserved charge, the simple result ${V}_{s}^{\mathcal{O}}(\lambda )={q}_{s}^{\mathrm{d}\mathrm{r}}(\lambda )$ holds true. Lastly, σ_s (λ) is a statistical factor determined by the nature of the excitations and it is different if quantum particles (bosons or fermions) or classical radiation or solitons are considered [112]. In our case, one has

$$\begin{equation}\sigma (\lambda )=\vartheta (\lambda )\enspace ,\qquad \varkappa > 0\enspace ;\qquad {\sigma }_{s}(\lambda )=1\enspace ,\qquad \varkappa < 0.\end{equation} \tag{ 66 }$$

Of course, the Eulerian correlation functions can also be accessed by first writing the quantum expression and then taking the semiclassical limit, consistently obtaining the same result. We now benchmark our analytical predictions against microscopic simulations. For numerical reasons clarified later on, rather than starting from a repulsive thermal state and drive the interactions in the attractive phase, we start by engineering suitable GGEs at the non-interacting point ϰ = 0 of the form

$$\begin{equation}{\left.\rho (\lambda )\right\vert }_{\varkappa ={0}^{+}}=\frac{1}{w\sqrt{\pi }}{\mathrm{e}}^{-{\lambda }^{2}/{w}^{2}}.\end{equation} \tag{ 67 }$$

With this convention, the total density of particles is fixed to unity and the parameter w is tuned to address different initial states. Notice that the root density is decaying fast for λ → ∞: this is not the case for thermal ensembles, which makes the discretization of the GHD equations more difficult. GGEs in the form equation (67) are easy to be implemented in the microscopic simulations, see appendix A. In figure 2 we investigate several protocols obtained by gently driving the system within the attractive phase ϰ < 0 starting from states (67). We consider linear ramps of the interaction ϰ(t) = −t/T_max and start by investigating the convergence to GHD as T_max is increased (figure 2(a)). We focus on the correlated density ${g}_{2}=\langle \vert \phi {\vert }^{4}\rangle /{(\langle \vert \phi {\vert }^{2}\rangle )}^{2}$ and observe a linear convergence in the ramp time 1/T_max, within the numerical error (see inset). Corrections to adiabaticity are of different nature: higher derivative terms in the GHD equation (62) are expected and a second order contribution $\sim {\partial }_{t}^{2}$ would lead to $\sim {T}_{\mathrm{max}}^{-1}\enspace$ corrections. Besides, finite-time protocols affect the bound state recombination happening at ϰ = 0⁻. Also in this case, one can envisage $\sim {T}_{\mathrm{max}}^{-1}\enspace$ corrections in view of the forthcoming argument. Following the quantum case discussed in reference [62], immediately after crossing ϰ = 0⁻, bound state recombination is allowed as long as the typical size of the bound state ∝|ϰ|⁻¹ exceeds the typical correlation length ζ_corr. In dimensionless units, this sets a recombination time ${t}_{\mathrm{B}\mathrm{S}}\propto {T}_{\mathrm{max}}\enspace {\zeta }_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}^{-1}$. Physically, only particles that remained close by during this time window are allowed to bind: let us consider two particles of rapidities λ and λ': on a timescale t_BS they drift apart by a distance t_BS(2λ − 2λ') (we neglect the interaction on this time scale and approximate the effective velocity with the bare velocity 2λ). Imposing the drift to be smaller than the bound state size at t_BS, i.e. ζ_corr, we obtained a coarse-grained resolution in rapidity space ${\Delta}\lambda \sim {\zeta }_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}^{2}/{T}_{\mathrm{max}}\enspace$. Thus, the Dirac deltas in rapidity space appearing in the maximum entropy ansatz will be broadened on a scale $\sim {\Delta}\lambda$, inducing overall $\propto {T}_{\mathrm{max}}^{-1}\enspace$ corrections. After having studied the approach to adiabaticity, we consider the protocol for different initial states focusing on g₂ (figure 2(b)) and on the density–density correlation function (65) (figure 2(c)), where we find excellent agreement. As expected, the gas is much more sensitive to attractive rather than repulsive interactions, as depicted in figure 3. There, we initialize the system at the non-interacting point and alternatively drive it within the attractive or repulsive phase. By retaining the same initial state, we explore the effects of the interactions. In particular, it is worthwhile to contrast the density correlation functions in the attractive phase with the repulsive phase (figures 3(b) and (c)). By increasing ϰ to negative values, the particles feel stronger attraction and convert part of their kinetic energy into interactions. Hence, a strong signal appears at small rays (i.e. small velocities). On the contrary, repulsive interactions accelerate the particles to larger velocities and reduce the peak of the correlation function. Finally, we address the original proposed protocol: in figure 4 we initialize the system in a thermal ensemble within the repulsive phase and adiabatically drive it to attractive interactions. Ideally, thermal states are simply a different initial state which can be captured by our method, but in practice the discretization of the GHD equations appears more challenging. The reason is in the shape of $\rho (\lambda ){\vert }_{\varkappa ={0}^{+}}$, which features polynomially decaying tails at large rapitidities $\sim {\lambda }^{-2}$, while it is very peaked at small rapitidies. As a consequence, a very large cutoff in the rapidity space must be used, while in the attractive phase the cutoff along the s direction greatly depends on λ. These features, in addition to the need of a fine discretization in rapidity space to correctly capture the integration kernels and the Eulerian shift of the GHD equation, worsen the convergence. We experienced that a λ-dependent discretization of the s space improves the convergence of the integral equations, but negatively affects the interpolation used in solving the GHD equation. Hence, we compromise with an independent, but inhomogeneous discretization of the λ and the s space: a judicious discretization choice is crucial to boost the convergence (see caption of figure 4). Overall, a much denser discretization in the λ space compared with the string space is needed for correctly describing the derivatives and integration kernels. In figure 4 we obtained a good agreement by considering different discretizations and then extrapolating (see caption for more details), while in figures 2 and 3 convergence was obtained without the need of any extrapolation and by employing a flat discretization in the λ and s space (roughly $\sim 800$ points have been used in the discretization).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this work, we constructed the integrable hydrodynamics of the one-dimensional attractive non linear Schrӧdinger equation. We then applied these findings to adiabatic interaction changes from the repulsive to the attractive phase, observing soliton production and engineering a strongly excited nonequilibrium phase of matter. The integrability of the NLS stabilizes the attractive phase out of equilibrium, which would be highly unstable in the presence of thermalization. We constructed the GHDs of the classical model by considering the proper semiclassical limit of the quantum Bose gas: the excellent agreement with ab initio microscopic simulations confirms the correctness of the procedure, but this strategy leaves open some interesting questions.

The first one concerns the absence of radiative modes, as one could have naively expected from the inverse scattering analysis. Let us stress that in statistical ensembles in the thermodynamic limit the distinction between solitons and radiations is blurred since one can accommodate for an arbitrary number (divergent in the thermodynamic limit) of zeroes in the imaginary plane of the scattering data (30). In this perspective, one could interpret radiation as a condensation of light solitons (small binding energy) on the real axis of the spectral parameter. Within the (λ, s) phase-space, this means that modes at s → 0 should be macroscopically occupied, which is not the case within our protocol. In fact, when ϰ = 0 is crossed over in an adiabatic manner, solitons of arbitrary size s are created with equal probability and s = 0 does not stand out as a special point. It is natural to wonder if radiative modes can be excited by different protocols.

A second question concerns the classical interpretation of the Planck constant h. Indeed, in contrast with the repulsive phase, we showed that the Planck constant does not completely disappear after having taken the semiclassical limit and plays an important role as a regulator in the entropy of the soliton gas. Interestingly, the h → 0 limit can be explicitly taken on the TBA equations obtained after entropy maximization (58), but not for the entropy itself. Indeed, by replacing the non-singular parameterization (59) in the semiclassical Yang–Yang entropy (51), at most $\sim \mathrm{log}\enspace s$ singularities (and thus are integrable) appear. As a consequence, the $\sim \mathrm{log}\enspace h$ term is not counterbalanced by a divergence in δ(h). Since it does not seem possible to explicitly remove the Planck constant from the semiclassical Yang–Yang entropy, the same term should be recovered from a pure classical construction of the thermodynamics. Retrospectively, we can give a phenomenological interpretation of the Planck constant as a infrared regularization of the classical theory. Indeed, NLS solitons have a width $\sim {(s\vert \varkappa \vert )}^{-1}$: when addressing the thermodynamics of a large, but finite system of length L, large solitons should be removed from the spectrum giving a cutoff s|ϰ| ≳ L⁻¹. In our results, this role is taken by the cutoff δ(h) ∝ h appearing in equation (51). However, notice that the validity of this crude argument is harmed by the need of taking the non-interacting limit and furthermore it does not explain the appearance of the $\sim \mathrm{log}\enspace h$ in the entropy. A complete understanding of these points passes through a complete classical analysis of the attractive phase. We leave this interesting question open for the future. Lastly, it would be interesting to extend our analysis to the SG model [27, 114–116]. This relativistic field theory is integrable both at the quantum and classical level and its spectrum is understood in terms of topological excitations (kinks) and bound states of them (breathers). The interest is twofold: first, by changing the interaction one forms or destroys kink's bound states, in strict analogy with the NLS case. Secondly, the thermodynamics of the classical SG, which in contrast to the NLS is stable also at equilibrium, has been thoroughly addressed in the literature also by means of semiclassical analysis [95, 99, 117]. On the other hand, the SG model and the NLS are closely connected, since the second can be obtained as the non-relativistic limit of the first [118–120]. However, to the best of our knowledge, an analogue $\sim \mathrm{log}\enspace h$ regulator does not appear in the SG literature. It would be interesting to further investigate this point.

We thank Žiga Krajnik, Enej Ilievski, Tomaž Prosen, Benjamin Doyon, Oleksandr Gamayun and Andrea De Luca for fruitful discussions. We are grateful to Žiga Krajnik, Enej Ilievski and Oleksandr Gamayun for useful comments on our manuscript. RK and JSC acknowledge support from the European Research Council (ERC) under ERC Advanced Grant 743032 DYNAMINT. AB acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy-EXC-2111-390814868.

The data that support the findings of this study are available upon reasonable request from the authors.

The microscopic simulations are performed by discretizing the NLS ϕ(x) → ϕ_j on N lattice sites and lattice space a, thus the system's length is L = aN. The limit of small lattice spacing is then numerically taken. In this respect, one discretize the NLS Hamiltonian as

$$\begin{equation}H\to a\sum\limits _{j}\left\{{a}^{-2}\vert {\phi }_{j+1}-{\phi }_{j}{\vert }^{2}+\varkappa \vert {\phi }_{j}{\vert }^{4}\right\},\end{equation} \tag{ A1 }$$

which results in the equation of motion

$$\begin{equation}i{\partial }_{t}{\phi }_{j}=-\frac{1}{{a}^{2}}({\phi }_{j+1}+{\phi }_{j-1}-2{\phi }_{j})+2\varkappa \vert {\phi }_{j}{\vert }^{2}{\phi }_{j}.\end{equation} \tag{ A2 }$$

Periodic boundary conditions are assumed. The evolution of the observables of interest is performed in two steps

• (a)  
  We randomly sample the initial condition ϕ_j (t = 0) from a known probability distribution. We consider two cases which are sampled in different ways:

  • 1.  
    For initial values ϰ > 0 we consider thermal initial states of inverse temperature β and chemical potential μ. Therefore, the target probability is $P[\left\{{\phi }_{j}\right\}]\propto \mathrm{exp}\left[-\beta (H-\mu a{\sum }_{j}\vert {\phi }_{j}{\vert }^{2})\right]$. Boltzmann distributions of this form are efficiently sampled by means of metropolis-hasting algorithms: basically, one builds a stochastic random walk on the field configuration $\left\{{\phi }_{j}\right\}\to \left\{{\phi }_{j}^{\prime }\right\}$, where the proposed update is accepted or rejected depending on $P[\left\{{\phi }_{j}^{\prime }\right\}]/P[\left\{{\phi }_{j}\right\}]$. This process is ergodic in the phase space and, after a sufficient number of initial steps, the average over the probability P can be replaced with an average over the metropolis time. These techniques are standard and the reader can refer to references [106–108] for further details.
  • 2.  
    By choosing ϰ = 0 as the initial point of the time evolution, we can conveniently initialize the state in arbitrary GGEs. Indeed, non-interacting GGEs are Gaussian states diagonal in the Fourier space

    $$\begin{equation}\langle {\tilde{\phi }}_{j}\rangle =0\qquad \langle {\tilde{\phi }}_{j}^{\ast }{\tilde{\phi }}_{{j}^{\prime }}\rangle =\frac{2\pi }{L}{\delta }_{j,{j}^{\prime }}\rho \left(\frac{2\pi }{L}{j}_{\mathrm{mod}(N,-N/2)}\right),\end{equation} \tag{ A3 }$$

    where ${\tilde{\phi }}_{j}$ is the field in the discrete Fourier space

    $$\begin{equation}{\phi }_{j}=\sum\limits _{{j}^{\prime }}\;{\mathrm{e}}^{i2\pi j{j}^{\prime }/N}{\phi }_{{j}^{\prime }}\end{equation} \tag{ A4 }$$

    and j_{mod(N,−N/2)} is the periodic repetition of j modulus N computed within the interval [−N/2, N/2]. Above, ρ(k) sets the initial root density in the limit ϰ → 0⁺ and can be set at will.
• (b)  
  Each random initial field configuration ϕ_j (t = 0) is then evolved according to the deterministic equation of motion (A2). In order to achieve stability, we trotterize the evolution on the time interval Δt in two steps. In the first step, we diagonally evolve in the real space according to

  $$\begin{equation}{\phi }_{j}(t+{\Delta}t/2)={\text{e}}^{-\text{i}{\Delta}t\varkappa (t)\vert {\phi }_{j}(t){\vert }^{2}}{\phi }_{j}(t).\end{equation} \tag{ A5 }$$

  In the second step, we instead evolve in the Fourier space

  $$\begin{equation}{\tilde{\phi }}_{j}(t+{\Delta}t)={\text{e}}^{-\text{i}\frac{2{\Delta}t}{{a}^{2}}[1-\mathrm{cos}(2\pi j)]}{\tilde{\phi }}_{j}(t+{\Delta}t/2).\end{equation} \tag{ A6 }$$

  This strategy explicitly enforces the particle number conservation. For a given lattice spacing a, we chose Δt by checking energy conservation of an initial random configuration in the case of a constant interaction ϰ. A posteriori, the convergence of the simulation is checked by considering different discretizations in space and time.

We start with discretizing the rapidity space and in case of attractive interactions the string space, additionally. We here choose a possibly inhomogeneous, but independent discretization in the rapidity and string spaces $[{\left\{{\lambda }_{i}\right\}}_{i=1}^{{N}_{\lambda }},{\left\{{s}_{\alpha }\right\}}_{\alpha =1}^{{N}_{s}}]$ which gives enough flexibility to compromise between numerical effectiveness and a sufficient resolution (further details on the concrete discretizaton can be found in the caption of figure 4). We moreover toss away points of the grid in the (λ, s)-plane where the filling fraction drops to zero. Within the attractive case, we use the non-singular reparameterization discussed in equation (59) and below.

Next, several integral equations need to be discretized for which we take the dressing (60) as an example in order to illustrate our proceedings. We use the grid in order to discretize equation (60) and get

$$\begin{equation}{s}_{\alpha }{t}_{{s}_{\alpha }}^{\mathbf{d}\mathbf{r}}({\lambda }_{i})={t}_{{s}_{\alpha }}({\lambda }_{i})-\sum\limits _{\beta =1}^{{N}_{s}}\enspace \sum\limits _{j=1}^{{N}_{\lambda }}\enspace {{\Phi}}_{{s}_{\alpha },{s}_{\beta }}({\lambda }_{i},{\lambda }_{j}){\bar{\vartheta }}_{{s}_{\beta }}({\lambda }_{j}){t}_{{s}_{\beta }}^{\mathbf{d}\mathbf{r}}({\lambda }_{j}).\end{equation} \tag{ B1 }$$

Here, ${{\Phi}}_{{s}_{\alpha },{s}_{\beta }}({\lambda }_{i},{\lambda }_{j})$ is the integration kernel analytically integrated over the domain surrounding (λ_j , s_β )

$$\begin{equation}{{\Phi}}_{{s}_{\alpha },{s}_{\beta }}({\lambda }_{i},{\lambda }_{j})=\underset{\frac{1}{2}({s}_{\beta }+{s}_{\beta -1}){}}{\overset{\frac{1}{2}({s}_{\beta }+{s}_{\beta +1})}{\int }}\mathrm{d}{s}^{\prime }\underset{\frac{1}{2}({\lambda }_{j}+{\lambda }_{j-1}){}}{\overset{\frac{1}{2}({\lambda }_{j}+{\lambda }_{j+1})}{\int }}\mathrm{d}{\lambda }^{\prime }{s}_{\alpha }^{-1}{\varphi }_{{s}_{\alpha }{s}^{\prime }}({\lambda }_{i}-{\lambda }^{\prime }).\end{equation} \tag{ B2 }$$

The integral above can be analytically performed. Notice that the integrand becomes singular for small interactions, pointing out the necessity of the discretization equation (B2) when compared with a more naive midpoint rule. A similar strategy is used in computing the force term and in the repulsive phase as well.

To start the time evolution determined by the GHD equations, we need to compute the initial state characterized by a root density. Within the repulsive phase, we initialize the system in a thermal state which we obtain upon solving the TBA equation (41) with a fixed chemical potential μ. Alternatively, we can construct a root density in form of a Gaussian (67) at the free point ϰ = 0⁺ which comes in handy to control the fat tails in rapidity space that appear in the thermal state. To enter the attractive phase, we use the root density at ϰ = 0⁺ as an input and maximize the entropy. We found that discretizing the integral equation (58) and then numerically solving them to yield exactly the input root density provides better accuracy rather than discretizing the analytical exact solution. The integral appearing in equation (58) is computed by a similar discretization strategy as the one used for the ϰ ≠ 0 dressing. In view of the constant cutoff in the string space employed in our discretization, we directly solve (58) by restricting the integral on a interval [0, S] and then discretize it on a finite grid. This strategy, compared with discretizing directly (52), allows to obtain an exact particle number conservation also on a finite discretization grid. The GHD equations are then solved in the filling space with the methods of characteristics. First, one observes that the GHD equation for the filling fraction (62) can be exactly recast as an implicit translation

$$\begin{equation}{\vartheta }_{s}(t+{\Delta}t,\lambda )={\vartheta }_{s}\left(t,\lambda -\underset{t{}}{\overset{t+{\Delta}t}{\int }}\mathrm{d}\tau {F}_{s}^{\mathrm{e}\mathrm{ff}}(\lambda ,\tau )\right),\end{equation} \tag{ B3 }$$

where, for simplicity, we consider the homogeneous case. The time integral of the force is computed with a midpoint rule, obtaining a second order algorithm in the time displacement. For details see references [46, 109].