# Observation of a generalized Gibbs ensemble in photonics

Alvise Bastianello, Alexey Tikan, François Copie, Stephane Randoux, Pierre Suret

► To cite this version:

Alvise Bastianello, Alexey Tikan, François Copie, Stephane Randoux, Pierre Suret. Observation of a generalized Gibbs ensemble in photonics. Physical Review A, 2026, 113 (1), pp.013514. 10.1103/xz2w-hndp. hal-05457000

HAL Id: hal-05457000

https://hal.science/hal-05457000v1

Submitted on 13 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Observation of a generalized Gibbs ensemble in photonics

Alvise Bastianello,¹,²,\* Alexey Tikan,³,\* Francois Copie,⁴ Stephane Randoux,⁴ and Pierre Suret⁴,†

¹Technical University of Munich, TUM School of Natural Sciences, Physics Department, 85748 Garching, Germany

²Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, 80799 München, Germany

³Laboratoire Temps-Fréquence, Université de Neuchâtel, Avenue de Bellevaux 51, Neuchâtel, Switzerland

⁴Univ. Lille, CNRS, UMR 8523 - PhLAM - Physique des Lasers Atomes et Molécules, F-59 000 Lille, France

In generic classical and quantum many-body systems, where typically energy and particle number are the only conserved quantities, stationary states are described by thermal equilibrium. In contrast, integrable systems showcase an infinite hierarchy of conserved quantities that inhibits conventional thermalization, forcing relaxation to a Generalized Gibbs Ensemble (GGE) – a concept first introduced in quantum many-body physics. In this study, we provide experimental evidence for the emergence of a GGE in a photonic system. By investigating partially coherent waves propagating in a normal dispersion optical fiber, governed by the one-dimensional defocusing nonlinear Schroedinger equation, we directly measure the density of states of the spectral parameter (rapidity) to confirm the time invariance of the full set of conserved charges. We also observe the relaxation of optical power statistics to the GGE’s theoretical prediction, obtained using the experimentally measured density of states. These complementary measurements unambiguously establish the formation of a GGE in our photonic platform, highlighting its potential as a powerful tool for probing many-body integrability and bridging classical and quantum integrable systems.

# I. INTRODUCTION

Strongly nonlinear systems and field theories out of equilibrium are among the most arduous challenges of modern physics. While determining every aspect of time evolution is a daunting problem, large scales see the emergence of universality and thermodynamics, dictating relaxation to the maximally entropic state compatible with conservation laws [1, 2]. When only the Hamiltonian is conserved, the system relaxes to a Gibbs thermal state. However, there are several known exceptions to this paradigm [3–5], such as integrable systems. Integrable systems are strongly interacting one-dimensional models with an infinite number of extensive conservation laws that constrain the dynamics and defy naive thermodynamics. Cutting-edge experiments [6–8] spurred two very active but seemingly disconnected communities in tremendously advancing understanding nonequilibrium, the first focusing on quantum systems [9, 10] and the other on classical partial differential equations (PDEs) [11]. Since emergent statistical properties are determined by symmetries and conservation laws, both communities independently rediscovered similar concepts from different angles.

The quantum community focused on quantum quenches [12], simple initial-value problems where a quantum state is unitarity evolved. Under minimal assumptions [13, 14], at late times an integrable system locally relaxes to a Generalized Gibbs Ensemble (GGE), first conjectured by Rigol et al. [15] in 2007. The GGE state $\hat{\rho} \propto e^{-\sum_n \mu_n \hat{Q}_n}$ generalizes conventional thermal Gibbs ensemble to the infinite number of conserved quantities, or charges, of integrable systems $\partial_t \hat{Q}_n = 0$. The

coefficients $\mu_n$ are determined by imposing charge conservation between the GGE and the initial conditions. The GGE concept opened to several groundbreaking theoretical achievements [16–21], and the GGE signature in correlation functions has been experimentally verified in cold atoms [22]. Experimental refinements further characterized GGEs, by recovering the root density or density of states [23]–equivalent to measuring all the conserved charges and thus uniquely identifying the GGE [24]–from atomic cloud expansions [25–30].

In parallel, Zakharov introduced the concept of integrable turbulence for PDEs in 2009, describing the evolution of random initial conditions in integrable systems such as the 1D nonlinear Schroedinger (NLS) or

<figure><img src="image_2.png" /><figcaption>FIG. 1. The GGE from integrable dynamics in optical fibers.— Panel (a): Experimental realization of the NLS through light propagation in nonlinear optical fibers. Panel (b): Evolution of a representative field configuration in the strongly nonlinear regime. Physical time and fiber lengths in the experiment are respectively proportional to space and time in the adimensional NLS equation (1).</figcaption></figure>

<footnote>\* These authors contributed equally to this work.</footnote>

<footnote>† pierre.suret@univ-lille.fr</footnote>
Korteweg–De Vries equations [31]. This seminal work has since inspired numerous theoretical studies [32–34] as well as experiments in optics [35, 36] and hydrodynamics [37–39]. Most experimental investigations have focused on statistical measurements, particularly the single-point probability density function (PDF) of the intensity [35, 36, 40, 41]. Numerical simulations and experiments have shown that, through integrable turbulence, the system eventually relaxes toward a statistically stationary state [31, 32, 36]. The natural assumption is relaxation to a GGE determined by the density of states in analogy with the quantum framework.

In this work, we present the experimental observation of relaxation to a GGE following a “quench” in a classical system, realized through the propagation of light in a single-mode optical fiber from random initial conditions as illustrated in Fig. 1(a). Light propagation in the laboratory frame is captured by the integrable defocusing NLS equation [42] $i2\partial_z A = \beta_2\partial_T^2 A - 2\gamma|A|^2 A$, where A is the slowly varying envelope of the optical field, whereas $\beta_2$ and $\gamma$ are the group velocity dispersion and Kerr nonlinearity of the fiber and z and T indicate the propagation distance along the fiber and the moving time frame. Through this work, we use adimensional units and write the NLS equation in the canonical form

$$i\partial_t\psi = -\partial_x^2\psi + 2|\psi|^2\psi, \quad (1)$$

where $\psi = A/\sqrt{P_0}$ and $P_0$ is the average optical power. The adimensional time t and position x are connected with the laboratory coordinates as $t = -\gamma P_0 z/2$, and $x = \sqrt{\gamma P_0/|\beta_2|}T$, see Appendix A for details. The GGE emerges as a statistical description of the field after a long time t equivalent to propagation through a long fiber, see Fig. 1(b).

Previous works experimentally [43] and theoretically [44] studied the intensity’s PDF obtained from initial configurations within a compact support, propagated after long times to distances much larger than its initial width. In contrast, the GGE emerges when the thermodynamic limit is taken first, and the long-time limit afterwards. We extract the density of states directly from the experimental data before and after evolution, demonstrating its conservation in time, and show it captures the emergent GGE by focusing on the PDF of the intensity. The intensity’s PDF $P(|\psi|^2)$ is the probability of measuring a certain value of $|\psi(x)|^2$ in a translationally invariant ensemble of fields: we compare the profiles evolved from the initial random field configurations with recent theoretical predictions on GGEs [45]. This measurement is challenging because probing the density of states requires an accurate single-shot recording of both phase and amplitude of the field. Our work places experimental platforms for classical integrable PDEs as a promising venue for investigations historically pertaining to the quantum realm, and poses a crucial milestone in establishing a fruitful dialogue between the communities working on quantum and classical many-body integrability, recently initiated with the parallelism [46] between

soliton gases [8, 47–49] and the more recent Generalized Hydrodynamics [50, 51], further strengthened by semi-classical limits of quantum integrability [45, 52–57].

## II. THE EXPERIMENTAL SETUP

The conceptual scheme of the experimental setup is illustrated in Fig. 1(a). The initial random conditions are prepared by shaping the spectrum of the amplified spontaneous emission (ASE) signal (Highwave Optical Technologies operating at central wavelength 1550 nm) with a programmable optical filter to a Gaussian profile. The optical power after a post-amplification stage has been adjusted with a variable attenuator. As a result, the initial condition corresponds to the superposition of independently and gaussianly distributed Fourier modes having random phase. This initial field distribution gives an exponential PDF of the optical power [32]. The random light is injected in 900 m-long single mode fiber with a group velocity dispersion coefficient $\beta_2 = 22 \text{ ps}^2/\text{km}$ and a nonlinear Kerr coefficient $\gamma = 3 \text{ W}^{-1}\text{km}^{-1}$. Our study requires the full field detection of non-periodic ultrafast optical signals, typically in the sub-picosecond regime. For this purpose, we employed a heterodyne time microscope (HTM) [35], which leverages space-time imaging analogies [58]. The HTM enables single-shot recording of both phase and power profiles of the optical signals with a temporal resolution of 250 fs over a temporal window of 40 ps. The details of the data acquisition process with HTM are provided in the Appendix A as well as in Refs. [35, 59].

The effective strength of the nonlinearity is tuned by controlling the adimensional momentum width $\Delta k = 2\pi\sqrt{\beta_2/(\gamma P_0)} \Delta\nu$, where $\Delta\nu$ is the optical bandwidth (in physical unit, Hz) of the initial signal. For large $\Delta k$, the kinetic (linear) term in the NLS equation (1) is dominant and the nonlinearity is weak. For small $\Delta k$, the nonlinear interaction gains importance. In this study, we used $\Delta\nu = 0.5 \text{ THz}$ and $P_0 = 3.5 \text{ W}$ for the weakly nonlinear case, and $\Delta\nu = 0.1 \text{ THz}$ and $P_0 = 3.6 \text{ W}$ for the strongly nonlinear case which corresponds to $\Delta k = 2.73$ and $\Delta k = 0.539$ for the normalized momentum, respectively, see Appendix A.

## III. THE DENSITY OF STATES

In translationally invariant and linear PDEs, normal modes are the Fourier components and their spectrum –which is conserved in time– is the density of states. In nonlinear, or interacting, integrable models normal modes become highly nonlinear functions of the field and the density of states is not the Fourier spectrum any longer. Instead of a simple Fourier transform, identifying normal modes requires the Inverse Scattering Transform (IST) [11] in classical PDEs, and the Bethe Ansatz [10] in quantum integrability. Here we briefly overview stan-
<figure><img src="image_3.png" /><figcaption>FIG. 2. Extracting the GGE from experimental data.— Panel (a): We show the profile of the norm and phase of one random field configuration after the weak and strong nonlinear evolution. To avoid boundary detection problems, we consider a central window (shadowed area $x \in [-8,8]$) and compute the density of states within it, shown in the rightmost panel in the respective colors. Panel (b): We show the conservation of the density of states by reporting its average values before and after the evolution, for weak and strong nonlinearities. When averaging, we consider $10^3$ experimental samples: the confidence interval (shaded areas) is obtained by considering the maximum spreading of partial averages done on 250 samples each.</figcaption></figure>

dard concepts of non-linear modes in integrable models and their density of states, and leave to the Appendices B and C a detailed and pedagogical overview in the quantum and classical cases respectively. Nonlinear modes in classical PDEs are generally radiation and solitons [11], whereas fermionic excitations emerge [23] in the quantum case. In both cases, modes are characterized by a spectral parameter $\lambda$ which in the quantum case, and through this work, is called the rapidity. These modes behave as ballistic quasiparticles showcasing only factorized elastic scattering [60], completely determined by the two-body scattering shift. In practice, their dynamics is reminiscent of a hard rods gas [61, 62], i.e. extended particles with elastic pairwise scattering, where the effective length depends on the rapidities of the colliding particles [63]. Each mode additively contributes to conserved charges: in the thermodynamic limit, one conveniently introduces a coarse grain mode density $\rho(\lambda)$ as $Q_n = L \int d\lambda q_n(\lambda)\rho(\lambda)$, where L is the large system size and $q_n(\lambda)$ the single-mode contribution to the charge. Since $\rho(\lambda)$ it parametrizes the charges' densities, which are thermodynamically-intensive quantities, it is self-averaging. In the classical and quantum cases, $\rho(\lambda)$ is commonly known as the density of states [48] or root density [23] respectively. Since $\rho(\lambda)$ determines all the charges, it is conserved in time and unambiguously fixes the GGE [24, 64]. In classical PDEs, the density of states is retrievable from the IST spectrum of field's configurations [65]: hereon, we consider the defocusing NLS. In this case, the spectrum is encompassed by a single radiative mode with real rapidity $\lambda$, whose density of states is computable as follows [45, 52, 56]. We focus on smooth field configurations in a large box of volume L centered in zero and introduce the auxiliary linear problem [11] $\partial_x\Psi(x) = U(x)\Psi(x)$, where $\Psi(x)$ is an auxiliary 2-component complex field. The $2 \times 2$ matrix $U(x) = \mathfrak{R}[\psi(x)]\sigma^x + \mathfrak{J}[\psi(x)]\sigma^y + \frac{\lambda}{2i}\sigma^z$ is the parallel-lem transporter and $\sigma^{x,y,z}$ are the canonical Pauli matrices. The auxiliary problem transforms the initial condi-

tion $\Psi(-L/2) = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ into $\Psi(L/2) = \begin{pmatrix} e^{-ix\lambda}a(\lambda) \\ b(\lambda) \end{pmatrix}$, whose first component $a(\lambda)$ generates the conservation laws $Q_n[11]$ and parametrizes the density of states [45, 52, 56] $\rho(\lambda) = \lim_{L\to\infty} \frac{1}{\pi L} \log|a(\lambda)|$, see also Appendix C. The density of states $\rho(\lambda)$ is then obtained tabulating the solution of the auxiliary problem for different values of $\lambda$. In the limit where the nonlinearity in Eq. (1) is negligible, the density of states approaches the Fourier spectrum and the rapidity becomes the Fourier wavevector, but they differ at strong nonlinearities. In Fig. 2, we show the density of states extracted from the experiment, both for the weakly nonlinear and strongly nonlinear regimes, and verify its conservation during time evolution. The density of states is extracted by discretizing and then numerically solving the auxiliary linear problem for different values of $\lambda$. See Appendix C for details on the numerical implementation and [66] for a commented code. Arbitrarily chosen examples of the recorded power and phase profiles after the propagation in the fiber are shown in Fig. 2a. The nonlinear evolution exhibits rich dynamics including dispersive shock waves [67] and dark soliton [68] formation, as shown in spatiotemporal diagram for similar conditions reported in Fig. 1(b). Interestingly, dark solitons do not explicitly contribute to the GGE, which is completely determined by the radiative mode only [45, 69]. However, kinetic equations for a dark solitons gas in the defocusing NLS have been derived in the framework of finite gap theory [70], raising interesting questions on the relation of the latter with the discussed GGE. The corresponding single-shot density of states are plotted in Fig. 2a (right panel) exhibiting a substantial shot-to-shot variation. To probe simultaneously conservation laws, we verify the time-independence of the density of states extracted from the experimental data. Field configurations at the initial and final time have been recorded for independent realizations, allowing us to verify the conservation of the density of states only on average. We analyzed 1000 data frames for both
the initial conditions and the signal after the propagation in the fiber, restricting the confidence interval to ≈26 ps window in the center of the frame, corresponding to $x \in [-8, 8]$ in rescaled units, to minimize the influence of reduced signal-to-noise ratio at the edges of the frame.

Importantly, fig. 2b demonstrates an excellent correspondence between the initial and final density of states $\rho(\lambda)$, confirming the conservation of charges and thus that we operate in the integrable regime, showing corrections to Eq. (1) coming from experimental imperfections (including losses) are negligible. The bandwidth of the initial Fourier modes determines the width of the density of states and the nonlinearity's strength. For weak interactions, the density of states is close to the momentum distribution and the two have similar width, whereas the two significantly depart in the strongly nonlinear regime.

## IV. EMERGENCE OF THE GGE

While the density of states remains unaffected through the dynamics, other observables showcase non-trivial evolution and eventual relaxation to a stationary state described by the GGE. For example, the PDF of the normalized power is initially an exponential $P_0(|\psi|^2) = e^{-|\psi|^2}$ as expected for the ensemble of independent Gaussian plane waves. At late times, under the effect of the nonlinear evolution, it relaxes to a new distribution determined by the GGE $P_{GGE}$. Our key result is showing the experimentally measured power's PDF relaxes to the corresponding theoretical prediction for GGEs derived in Ref. [45], using semiclassical limits [52] of previous results in quantum integrability [71, 72]. The PDF of the intensity on a GGE has the integral representation

$$P_{GGE}(|\psi|^2) = 16 \int_0^\infty d\theta J_0(8\sqrt{\theta}|\psi|^2) e^{-\frac{16}{\pi} \int_0^{\sqrt{\theta}} d\tau \tau S(\tau)}, \quad (2)$$

where $J_0(x)$ is a modified Bessel function of the first kind, $S(\tau) \equiv \int d\lambda s_\tau(\lambda)$, and $s_\tau(\lambda)$ solves the auxiliary integral equation

$$[\vartheta^{-1}(\lambda) + i0^{+}]s_{\tau}(\lambda) = 1 + \int \frac{d\lambda'}{2\pi} \frac{2}{\lambda - \lambda'} (4i\tau - \partial_{\lambda'})s_{\tau}(\lambda'). \quad (3)$$

Above, the singular integral is meant within a principal value regularization, and we introduce other standard quantities from quantum integrability [23] as the filling fraction $\vartheta(\lambda) = \rho(\lambda)/\rho^t(\lambda)$, and the total root density $\rho^t(\lambda) = \frac{1}{2\pi} - \int \frac{d\lambda'}{2\pi} \frac{2}{\lambda - \lambda'} \partial_{\lambda'} \rho(\lambda')$. The regularization $i0^+$ is needed for Eq. (3) to be well-defined, and it results in zeroes in $e^{-\frac{16}{\pi} \int_0^{\sqrt{\theta}} d\tau S(\tau)}$: details on the practical computations of the PDF are provided in Appendix D, see also Zenodo [66] for a working code. As we show in the Appendix D, Eq. (2) recovers the expected exponential distribution for weak non-linearities, but it sensibly departs from it at stronger non-linearities. In Fig. 3, we compare experimental measurements of the PDF with theoretical

<figure><img src="image_4.png" /><figcaption>FIG. 3. Relaxation of the intensity PDF to the GGE.— We show the PDF of the intensity field $|\psi|^2$ extracted from the initial data on the initial conditions and after evolution (red line and blue line respectively), in linear (a,b) and log (c,d) scale. We compare PDF data with the GGE result (black line) (2) obtained from the density of states extracted from experimental data, see Fig. 2. We compare with the exponential decay from Gaussian independent plane waves as reference (dashed line). As before, we show the average over $10^3$ independent samples, the confidence interval is the maximum spreading from partial averages of 250 samples each. The experimental PDF is extracted from the same window of data $x \in [-8, 8]$ used for the density of states in Fig. 2, and the bin size is 0.03.</figcaption></figure>

predictions, having the experimental density of states values (Fig 2b) as the input. The PDF was computed on the same 1000 independent field configurations previously used in Fig 2b, and considering spatial averaging within the same central window. The confidence interval is estimated from the maximum spreading of partial averages on four groups of 250 samples each. The theory uncertainty is analogously estimated: we extract the density of states for each group of 250 samples, compute the relative PDF and consider the maximum spreading. The PDF solid line is obtained from the density of states average over all field configurations. Fig. 3a depicts the weakly nonlinear case, where the PDF of the stationary state only slightly deviates from the initial exponential distribution (red line). The discrepancy between the experimental data and the prediction at low power levels is mainly attributed to the finite resolution of the HTM. Fig. 3b illustrates the strongly nonlinear case, where the PDF undergoes a substantial transformation due to nonlinear propagation, with a non-monotonous behavior. Also in this regime, experimental data agree well with GGE theoretical predictions. Altogether, the combined measurements provided in Fig. 2 and Fig. 3 show the emergence of a GGE stationary state, determined by the initial conditions through the conservation laws.
V. DISCUSSION

Photonic devices, particularly optical fibers, serve as precise and versatile platforms for studying nonlinear waves. Our results demonstrate that photonics can address fundamental questions traditionally explored in quantum many-body physics. Specifically, we provide a full characterization of the classical analog of a “quantum quench” [12] within the integrable classical defocusing NLS equation. We investigated the emergence of a GGE from a random-wave initial condition by measuring the density of states, analogous to the rapidity distribution in quantum systems, and confirming the relaxation of the experimental PDF of $|\psi|^2$ to the theoretical GGE prediction [45]. The characterization of the dynamics of the quantum counterpart of the NLS equation, the Lieb-Liniger model [73, 74], attracted considerable experimental effort in cold atoms [75]. The density of states, or rapidity distribution as it is often denoted in these experiments, is indirectly measured through atomic cloud expansions [25–30], and fluctuations of the number of particles [76, 77], i.e. the quantum analog of the PDF of $|\psi|^2$, can be measured with a finite-interval resolution. In contrast, photonic platforms enable direct access to the density of states and high-resolution power fluctuation measurements via local detection of the complex field. These high-precision measurements of the field enable a complete characterization of the state, whereas ultracold gases have access only to limited information: depending on the experimental platform, either coarse grain imaging of $|\psi(x)|^2$ or momentum distributions. Importantly, the destructive measurement process in quantum systems prevents the observation of two-time correlation functions, which can instead be accessed in photonics and have a key-role in the study of integrable systems [78, 79].

Notably, theoretical predictions for fluctuations of $|\psi|^2$ [45] were not obtained directly in the NLS, but through semiclassical limits of results in Lieb-Liniger [71, 72]. The study of semiclassical limits in quantum integrability has recently gained considerable momentum [45, 52–57] in addressing unresolved classical problems. Our work elevates this quantum-classical correspondence from theory to experiments, establishing photonics as a powerful tool to complement and potentially surpass cold atom systems in probing the statistical mechanics and hydrodynamics of integrable models.

ACKNOWLEDGEMENTS

AB acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy–EXC–2111–390814868. AB is grateful to G. Del Vecchio Del Vecchio, A. De Luca and G. Mussardo for collaboration in Ref. [45]. We are grateful to Jerome Dubail for useful comments on the manuscript. FC, SR and PS acknowledge support by the Agence

Nationale de la Recherche through the SOGOOD (ANR-21-CE30-0061) project, the LABEX CEMPI project (ANR-11-LABX-0007), the Ministry of Higher Education and Research, Hauts de France council and European Regional Development Fund (ERDF) through the Nord-Pas de Calais Regional Research Council and the European Regional Development Fund (ERDF) through the Contrat de Projets Etat-Région (CPER Photonics for Society P4S). FC, SR and PS thank the Centre d'Etudes et de Recherche Lasers et Application (CERLA) for technical support.

DATA AND INFORMATIONS AVAILABILITY.

Data, data analysis, and simulation codes are available upon reasonable request on Zenodo [66].

Appendix A: Experimental data analysis

**Time Microscope signal description.** The phase and power profiles are reconstructed from the 2-D snapshot (see in Fig. 4) recorded with the Heterodyne Time Microscope [35, 59]. The vertical line $y$ at position $x$ of a 2-D snapshot detected with the sCMOS camera of HTM reflects the interference between the non-collinear signal and reference beams resulting in the periodic pattern at a given time. The phase of the reference is considered to be a constant during the snapshot time window. Thus, the relative position of the interference patterns contain the information about the phase of the signal. The intensity distribution in the 2D snapshot can be formally expressed as:

$$I(x, y) = I_r(x, y) + I_s(x, y) + 2\sqrt{I_r(x, y)I_s(x, y)} \cos[k_y y + \phi(x)], \quad (A1)$$

where $I_r(x, y)$ transverse reference profile, and $I_s(x, y)$ - the transverse the signal profile, $2\pi/k_y$ is the spatial period of the interference pattern along $y$, observed by the sCMOS camera, $\phi(x)$ corresponds to the phase of the complex field envelope. The angle between the reference and signal beams is chosen to ensure sufficient fringe density, satisfying:

$$\int \sqrt{I_r(x, y) \cdot I_s(x, y)} \cos[k_y y + \phi(x)] dy \approx 0. \quad (A2)$$

**Optical power profile extraction.** For each frame, the optical power profile $P(x)$ of the signal is computed as:

$$P(x) = \int I(x, y) dy - \int I_r(x, y) dy, \quad (A3)$$

where $P(x)$ is expressed in arbitrary units. The reference intensity is detected in the absence of the signal beam
<figure><img src="image_5.png" /><figcaption>FIG. 4. An example of the experimentally recorded heterodyne time microscope frame and corresponding power and phase profile extraction. (Top) A representative 2-D snapshot recorded experimentally with the heterodyne time microscope, captured by the sCMOS camera. The spatial interference fringes arise from the non-collinear interaction of signal and reference beams at the detection plane, encoding both amplitude (brightness) and phase (relative fringe position) information in the spatial domain. (Middle) Corresponding temporal power profile extracted from the snapshot, obtained by integrating the intensity distribution along the vertical coordinate. (Bottom) Reconstructed temporal phase profile derived through spatial Fourier analysis of the recorded fringe pattern, which contains the relative phase information of the signal beam, with the reference beam's phase assumed constant across the measurement window.</figcaption></figure>

and averaged over $5 \times 10^4$ frames to exclude the fluctuations. The power is normalized as $P(x)/\langle P(x) \rangle$ or scaled to $[P(x)/\langle P(x) \rangle] \times P_0$, where $P_0$ is the average power launched inside the fiber. $\langle P(x) \rangle$ is computed with $5 \times 10^3$ frames to exclude the influence of injected power fluctuation during the experiment. We note that since the average power envelope $\langle P(x) \rangle$ has a limited duration, we observe the degradation of the signal-to-noise ratio at the edge of the extracted signal power $P(x)/\langle P(x) \rangle$. For this reason, for comparison with the theoretical predictions, we used a window reduced in both $x$ (typically to 200 pixels) and $y$ (100 pixels) directions. Our time calibration measurements gave one pixel in $x$: 79 fs correspondence.

**Phase profile extraction.** The relative phase $\phi(x)$ between the signal and reference beams is deduced from the interference fringes. Taking the Fourier transform of the frame in the vertical direction, we extract the phase of the Fourier component that corresponds to the fringes period. Thus the expression for the phase extraction takes the following form:

$$ \phi(x) = \arg(\max(\text{FFT}[I(x, y)]_y)), \quad (A4) $$

where the maximum amplitude is taken avoiding the low

frequency Fourier components.

**Computing the probability density function.** To compute the PDF, we analyzed 1000 frames. The confidence interval corresponds to partial averages on four groups equal groups of 250 frames. The data window is restricted to 26 ps in the center of the frame. We note that the experimental study of the initial conditions and signal after the fiber propagation are recorded separately, in some cases on different days. We, however, ensured the similar performance of the HTM by providing the calibration procedures described in [35].

**Normalization of the NLS.** We start with the dimensional form of the NLSE widely used in fiber optics [42]:

$$ \frac{\partial A}{\partial z_{\text{lab}}} = -i \frac{\beta_2}{2} \frac{\partial^2 A}{\partial t_{\text{lab}}^2} + i \gamma |A|^2 A. \quad (A5) $$

Where $z_{\text{lab}}$ and $t_{\text{lab}}$ are space and time coordinates in the laboratory. The normalization of the NLS used in the main text is obtained by defining

$$ \psi = \frac{A}{\sqrt{P_0}} \qquad t = -\frac{1}{2}\gamma P_0 z_{\text{lab}}; \qquad x = \sqrt{\frac{\gamma P_0}{\beta_2}} t_{\text{lab}}, \quad (A6) $$

where $\gamma$ is the nonlinear Kerr coefficient, $\beta_2$ is the group velocity dispersion and $P_0$ - the average optical power. In these units, the average intensity is normalized to unity $\langle |\psi|^2 \rangle = 1$, and the field evolves with the NLS equation

$$ i\partial_t \psi = -\partial_x^2 \psi + 2|\psi|^2 \psi. \quad (A7) $$

The initial random plane wave distribution is described by a gaussian ensemble for the field, most conveniently expressed in the Fourier space $\tilde{\psi}(k) = \int dx e^{ikx} \psi(x)$. The Fourier components are gaussianally distributed with zero mean $\langle \tilde{\psi}(k) \rangle = 0$ and variance $\langle \tilde{\psi}^*(k) \tilde{\psi}(q) \rangle \propto \delta(k-q)e^{-k^2/\Delta k^2}$, the prefactor is determined by the normalization condition $\langle |\psi|^2 \rangle = 1$ and the width of the distribution is

$$ \Delta k = 2\pi \sqrt{\frac{\beta_2}{\gamma P_0}} \Delta \nu, \quad (A8) $$

where $\Delta \nu$ is the experimental signal bandwidth.

## Appendix B: Overview of quantum integrability and thermodynamic Bethe ansatz

In this Appendix, for the sake of completeness, we provide a brief overview on those aspects of the quantum Lieb-Liniger model that are relevant for our scope. A more extensive and in-depth discussion is left to the literature mentioned through the Appendix.

**The microscopic solution: coordinate Bethe ansatz and Bethe equations.**— The determination
of the rapidity eigenstates $|\{\lambda_i\}_{i=1}^N\rangle$ is most conveniently framed within the coordinate Bethe ansatz. For a in-depth and pedagogical discussion of this Appendix, see Ref. [80]: we consider the Hamiltonian in first quantization within the $N$ particles sector $\hat{H} = \sum_i -\partial_{x_i}^2 + c \sum_{i \neq j} \delta(x_i - x_j)$, where $c$ parametrizes the interactions, and look for eigenstates. We consider the bosonic trial wavefunction $\phi_{\{\lambda_i\}_{i=1}^N}(x_1, \dots, x_N)$ and the ansatz

$$ \phi_{\{\lambda_i\}_{i=1}^N}(x_1 < \dots < x_N) \propto \sum_{\mathcal{P}} A(\mathcal{P}) \prod_{j=1}^N e^{i\lambda_{\mathcal{P}_j} x_j}, \quad (B1) $$

where the symmetric continuation over all the possible reordering of the coordinates $x_j$ is assumed. The summation is over all the permutations $\mathcal{P}$ of $N$ indexes $\{1, \dots, N\}$ and $A(\mathcal{P})$ are coefficients to be determined. When the coordinates $x_j$ are all different, the interaction term in the Hamiltonian is ineffective and Eq. (B1) reduces to a combination of plane waves, which are eigenstate of the derivative part of the Hamiltonian. Requiring $\phi_{\{\lambda_i\}_{i=1}^N}$ to be an eigenstate also at coincident coordinates leads to constraints on the coefficients $A(\Pi_{j,j+1}\mathcal{P}) = S(\lambda_{\mathcal{P}_j} - \lambda_{\mathcal{P}_{j+1}})A(\mathcal{P})$, with $\Pi_{j,j+1}$ the permutation of the indexes $j$ and $j+1$, and $S(\lambda) = \frac{\lambda+ic}{\lambda-ic}$ the so-called scattering matrix: in the LL model, $S(\lambda)$ is just a complex number with modulus one. An explicit solution for $A(\mathcal{P})$ can be found, but it is not crucial for our scope, and $\phi_{\{\lambda_i\}_{i=1}^N}$ (B1) becomes a true eigenvector of the Hamiltonian and of all the extensive conserved charges $\hat{Q}_n|\{\lambda_i\}_{i=1}^N\rangle = (\sum_{j=1}^N q_n(\lambda_j))|\{\lambda_i\}_{i=1}^N\rangle$. The charge eigenvalue $q(\lambda)$ associated to the Hamiltonian is $\lambda^2$, while for the momentum and number of particles one has $\lambda$ and 1 respectively. In general, charge eigenvalues of the form $q_n(\lambda) = \lambda^n$ are associated to the extensive conservation laws.

When considering thermodynamics, it is crucial having a finite density of particles. We assume the quantum particles live on a interval of length $L$ and ask for periodic boundary conditions (PBC) for the wavefunction: other boundary conditions, for example open boundary conditions, lead to the same thermodynamics when $L \to \infty$. Imposing PBC on Eq. (B1) leads to the Bethe equations [81], a set of nonlinear equations that quantize the rapidities

$$ e^{iL\lambda_j} = \prod_{k \neq j} S(\lambda_j - \lambda_k), \quad j, k \in \{1, \dots, N\}. \quad (B2) $$

All the solutions to the Bethe equations are associated to physical eigenstates Eq. (B1), provided the rapidities are chosen all different $\lambda_j \neq \lambda_k$: in case of equal rapidities, the wavefunction vanishes and so those configurations may be excluded. In the thermodynamic limit, complex solutions to the Bethe equations are associated with bound states and their presence depends on the sign of the interactions: for the repulsive case $c > 0$ we are interested in, only real solutions are important.

## GGEs and rudiments of thermodynamic Bethe

ansatz.— Armed with the Bethe equations, one can forget about the explicit form of the eigenstates (B1) in building the thermodynamics within the thermodynamic Bethe ansatz (TBA) [23], which we briefly overview (see also Ref. [80]). First, one considers the Bethe equations in logarithmic form and parametrizes the rapidities in terms of Bethe integers $I_j \in \mathbb{Z}$ as

$$ \frac{I_j}{L} = \frac{\lambda_j}{2\pi} - \frac{1}{2\pi L} \sum_{k \neq j} \Theta_q(\lambda_j - \lambda_k), \quad (B3) $$

where one defines the scattering phase as $\Theta_q(\lambda) = -i \log S(\lambda)$: the eigenstates $|\{\lambda_i\}_{i=1}^N\rangle$ are now in one-to-one correspondence with all the set of distinct integers $\{I_j\}_{j=1}^N$. Let us now consider the GGE partition function $\mathcal{Z} = \text{Tr}[e^{-\sum_n \beta_n \hat{Q}_n}]$, which one rewrites as an explicit sum over the eigenstates parametrized either with the rapidities, or the Bethe integers. The second choice is more convenient for further manipulations $\mathcal{Z} = \sum_N \sum_{\{I_j\}_{j=1}^N} e^{-\sum_{j=1}^N \sum_n \beta_n q_n(\lambda_j)}$, where $I_j$ and $\lambda_j$ are connected through the Bethe equations. We are now ready to move to the thermodynamic limit and, upon defining $y_j = I_j/L$, introduce a counting function $\sigma(y)$ which counts the number of $y_j$ within an interval of size $dy$, i.e. $Ldy\sigma(y) = \{\#y_j \in [y - \frac{dy}{2}, y + \frac{dy}{2}]\}$. The sum in the definition of the partition function is now readily converted into a path integral over the possible choices of the coarse grain counting function $\sigma$. When doing so, one should introduce a proper entropic weight due to the possible microscopic rearrangement of $y_j$ that do not affect $\sigma(y)$. Since $y_j$ behaves as hard-core particles, this eventually results in a Fermi-Dirac statistics [80]

$$ \mathcal{Z} \approx \int \mathcal{D}\sigma \exp \left[ L \int dy f(\sigma(y)) - L \int d\lambda \sum_n \beta_n q(\lambda(y)) \sigma(y) \right] \quad (B4) $$

where $f(x) = -x \log x - (1-x) \log(1-x)$ and one introduces the $\lambda(y)$ function from the coarse grain limit of the logarithmic Bethe equations $y = \frac{\lambda(y)}{2\pi} - \frac{1}{2\pi} \int dy' \Theta_q(\lambda(y) - \lambda(y'))\sigma(y')$. As a last step, in Eq. (B4) it is convenient to change variable $y \to \lambda$ by introducing the total root density $\rho^t(\lambda)$ defined through the Jacobian of the transformation $\rho_q^t(\lambda) \equiv \frac{dy}{d\lambda}$, and the root density $\rho_q(\lambda) = \rho_q^t(\lambda)\sigma(y(\lambda))$. After this change of variables, the GGE partition function is rewritten as

$$ \mathcal{Z} \approx \int \mathcal{D}\rho \exp \left[ L S[\rho_q] - L \int d\lambda \left( \sum_n q_n(\lambda) \right) \rho_q(\lambda) \right], \quad (B5) $$

where we define the entropy density

$$ S[\rho_q] = \int d\lambda \rho_q^t f(\vartheta_q(\lambda)), \quad (B6) $$

and the filling fraction $\vartheta_q(\lambda) \equiv \rho_q(\lambda)/\rho_q^t(\lambda)$. Computing the path integral (B5) is a daunting task, but in the thermodynamic limit $L \to +\infty$ only the saddle point matters.
Looking for the saddle point finally leads to integral equations connecting the Lagrange multipliers and the filling function (and thus the root density)

$$f'(\vartheta(\lambda)) = \sum_n \beta_n q_n(\lambda) + \int \frac{d\lambda'}{2\pi} \varphi_q(\lambda - \lambda') [\vartheta_q(\lambda)f'(\vartheta_q(\lambda)) - f(\vartheta_q(\lambda))]. \quad (B7)$$

Above, the scattering kernel is defined as $\varphi_q(\lambda) \equiv \partial_\lambda \Theta_q(\lambda)$, and in Lieb Liniger it reads $\varphi_q(\lambda) = -\frac{2c}{c^2+\lambda^2}$ [73].

**The density moments.**— The TBA equation (B7) is a cornerstone in describing thermodynamics and more in general GGEs, and it already allows for a direct computation of certain observables. For example, apart from more naive expectation values of local charges $\langle \hat{Q}_n \rangle = L \int d\lambda q_n(\lambda) \rho(\lambda)$, their connected correlators can be computed through derivatives of the partition function $\langle \hat{Q}_n \hat{Q}_{n'} \rangle_{\text{conn.}} = -\frac{\partial^2}{\partial \beta_n \partial \beta_{n'}} \log \mathcal{Z}$. In principle, the knowledge of $\rho(\lambda)$ and, through (B7), of the whole GGE determines the expectation values of all local observables. However, obtaining analytical formulas is a formidable task. Recently, close integral expressions for the density moments $\langle (\hat{\psi}^\dagger)^n (\hat{\psi}^\dagger)^n \rangle$ have been computed in Ref. [71, 72], and are encoded in the following generating function

$$1 + \sum_{n=1}^{\infty} Y^n \frac{(2c)^n}{n!} \langle (\hat{\psi}^\dagger)^n (\hat{\psi}^\dagger)^n \rangle = \exp \left( \frac{1}{\pi} \sum_{n=1}^{\infty} Y^n \mathcal{G}_n \right), \quad (B8)$$

from which the moments can be obtained upon expanding both sides in powers of the auxiliary variable $Y$. Above, $\mathcal{G}_n = \frac{c^{2n-1}}{n} \int d\lambda \vartheta_q(\lambda) \xi_{2n-1}^{\text{dr}}(\lambda)$, where the dressing operation $\xi_j(\lambda) \to \xi_j^{\text{dr}}(\lambda)$ is in general defined as the solution of the following integral equation

$$\xi_j^{\text{dr}}(\lambda) = \xi_j(\lambda) - \int \frac{d\lambda'}{2\pi} \varphi_q(\lambda - \lambda') \vartheta_q(\lambda') \xi_j^{\text{dr}}(\lambda'), \quad (B9)$$

and the auxiliary functions $\xi_j(\lambda)$ are determined by recursive integral equations

$$\begin{aligned} \xi_{2n}(\lambda) &= \int \frac{d\lambda}{2\pi} \vartheta_q(\lambda) \{\Gamma(\lambda - \lambda') [2\xi_{2n-1}^{\text{dr}}(\lambda') - \xi_{2n-3}^{\text{dr}}(\lambda')] + \\ & \qquad \varphi_q(\lambda - \lambda') \xi_{2n-2}^{\text{dr}}(\lambda')\}, \\ \xi_{2n+1}(\lambda) &= \delta_{n,0} + \int \frac{d\lambda'}{2\pi} \vartheta_q(\lambda') \vartheta_q(\lambda') \{\Gamma(\lambda - \lambda') \xi_{2n}^{\text{dr}}(\lambda') + \\ & \qquad \varphi(\lambda - \lambda') \xi_{2n-1}^{\text{dr}}(\lambda')\}, \end{aligned} \quad (B10)$$

with $\Gamma(\lambda) = -\frac{\lambda}{c} \varphi(\lambda)$ and $\xi_{n \le 0}(\lambda) = 0$. These recursive equations can be efficiently solved numerically [72]. This result on the density moments on a GGE in the quantum model gives access to the classical PDF through semiclassical limits, see Appendix D.

## Appendix C: Overview of classical integrability and thermodynamic Bethe ansatz

In this Appendix, we provide a short overview of the classical inverse scattering for the defocusing NLS and its connection with thermodynamics.

**Lax pairs and monodromy matrix.**— The NLS equation emerges as the compatibility condition for an auxiliary linear problem [11]. Let us define the $\lambda$-dependent Lax pair

$$U_{\lambda}(t, x) = (\psi^{*} \sigma^{+} + \psi \sigma^{-}) + \frac{\lambda}{2i} \sigma^{z}, \quad (C1)$$

$$V_{\lambda}(t, x) = \left( i|\psi|^2 + \frac{i\lambda^2}{2} \right) \sigma^z - [\lambda\sigma^x + i\partial_x\psi^*\sigma^+ - i\partial_x\psi\sigma^-],$$

with $\sigma^{\pm} = (\sigma^x \pm \sigma^y)/2$, and introduce the auxiliary field $\Psi(t, x) \in \mathbb{C}^2$ obeying the following linear differential equation

$$\partial_x \Psi(t, x) = U_{\lambda}(t, x) \Psi(t, x), \quad \partial_t \Psi(t, x) = V_{\lambda}(t, x) \Psi(t, x). \quad (C2)$$

The consistency of the auxiliary problem (C2) requires the zero-curvature condition $\partial_\lambda U_\lambda - \partial_x V_\lambda + [U_\lambda, V_\lambda] = 0$ which is satisfied if and only if the field $\psi$ obeys the NLS equation. Describing the NLS as a compatibility condition for a Lax pair manifests integrability and the presence of conservation laws. Below, we need to discuss separately two different boundary conditions for the NLS field: (i) the “whole line problem” where the field $\psi(t, x)$ vanishes at $|x| \to \infty$ faster than any power law, and (ii) the periodic problem $\psi(t, x) = \psi(t, x + L)$. A thermodynamic description of the state can be achieved combining these two problems.

**The whole-line problem.**— Let us define the spatial propagator for the auxiliary problem (C2) $T(t, x, y) = P\exp[\int_x^y dz U_\lambda(t, z)]$, which solves the simple 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)$. It is immediate to check that, whenever the interval $(x, y)$ is chosen sufficiently large in such a way the interval contains the whole non-trivial part of the field $\psi$, then $\partial_t \text{Tr}[T_\lambda(t, x, y)] = 0$ for any $\lambda$. Since the field $\psi$ initially vanishing outside of the interval $(x', y')$ can propagate beyond the initial interval during time evolution, we now send the extrema to infinity. However, for practical purposes it should be remembered that the considerations below apply for a finite interval as well, as long as the field vanishes out of it. We then define the transfer matrix on the whole line as [11, 52]

$$\mathcal{T}_{\lambda}(t) = \lim_{x \to \infty} e^{\frac{-\lambda x}{2i}} \sigma^z T_{\lambda}(t, -x, x) e^{-\frac{\lambda x}{2i}} \sigma^z, \quad (C3)$$

where the additional oscillating terms counterbalance the asymptotic behavior of $U_\lambda(t, x)$. Based on the symmetries of the linear problem only, the transfer matrix can be parametrized as

$$\mathcal{T}_{\lambda}(t) = \begin{pmatrix} a(\lambda) & b^{*}(\lambda) \\ b(\lambda) & a^{*}(\lambda) \end{pmatrix}, \quad (C4)$$
where for the sake of notation, the time-dependence of the matrix elements is left implicit. We have already seen as $\partial_t \text{Tr}[\mathcal{T}]_\lambda(t) = \partial_t [a(\lambda) + a^*(\lambda)] = 0$, but for vanishing boundary conditions it holds the stronger statement $\partial_t a(\lambda) = 0$, as it can be immediately checked from the equation for the time evolution of the propagator $T_\lambda(t, x, y)$. The coefficient $a(\lambda)$ is conserved for every $\lambda$, and it can be explicitly shown [11] the Laurent expansion of its logarithm generates the extensive charges $\log a(\lambda) \stackrel{\lambda \to \infty}{=} \sum_n \frac{1}{\lambda^{n+1}} Q_n$. We will now use this fact and the analiticity properties of $a(\lambda)$ to recover the classical counterpart of the quantum root density: an in depth discussion is left to the literature [11, 52], here we summarize the main facts. The coefficient $a(\lambda)$ can be analytically continued in the complex plane and $\lim_{|\lambda| \to \infty} a(\lambda) \to 1$: in the focusing regime, $a(\lambda)$ can have zeroes and poles in the complex plane, which are associated with solitons, while these are absent in the defocusing regime of interest for us. Therefore, one can use the Kramers-Kronig relations and parametrize $a(\lambda)$ as

$$a(\lambda) = \exp \left[ \int \frac{d\lambda'}{\pi i} \frac{\log |a(\lambda')|}{\lambda' - \lambda - i0^+} \right]. \quad (C5)$$

By taking the logarithm, and considering the Laurent expansion one reaches the identification $Q_n = \int d\lambda \lambda^n \frac{\log|a(\lambda)|}{\pi}$. A final identification of $\log|a(\lambda)|$ with the classical root density passes through the discussion of the periodic problem.

**The periodic problem.**— We here assume that the NLS field is periodic $\psi(t, x + L) = \psi(t, x)$, for further details see [52, 82]. As the auxiliary problem is linear and thus symmetric under rescaling of $\Psi(t, x)$, one does not ask the auxiliary field to be periodic, but allows for more freedom allowing periodicity modulus an overall complex constant $\Psi(t, x) = A\Psi(t, x + L)$ with unitary norm $|A| = 1$. A different choice $|A| \neq 1$ would lead to not normalizable solutions for $\Psi$ in the thermodynamic limit, and are thus excluded. As we see below, $A$ is further constrained. We introduce the PBC transfer matrix from the propagator as $\mathcal{T}_\lambda^{\text{PBC}} = T_\lambda(t, 0, L)$: notice that, in contrast with Eq. (C3), further oscillating terms are not included. In the periodic problem, the transfer matrix can still be parametrized as in Eq. (C4), and therefore $[\mathcal{T}_\lambda]_{11} = [\mathcal{T}_\lambda]_{22}^*$: this leaves as only possibility for the boundary conditions $A = \pm 1$. Notice that while the trace of $\mathcal{T}_\lambda^{\text{PBC}}$ is still a conserved quantity, the diagonal entries are not separately conserved. The periodic boundary conditions $\mathcal{T}_\lambda^{\text{PBC}} = \pm \text{Id}$ can be fulfilled only for discrete values of the spectral parameter $\lambda$, and leads to the classical counterpart of the Bethe equations introduced in Appendix B.

**The emergence of thermodynamics.**— We can now combine the knowledge of the whole-line and periodic problem to take the thermodynamic limit, see Ref. [52] for an extensive discussion. We consider a fundamental

interval $\mathcal{I}_1 = [-L/2, L/2]$, and a smaller one $\mathcal{I}_2 \subset \mathcal{I}_1$ defined as $\mathcal{I}_2 = [-L/2 + \ell/2, L/2 - \ell/2]$. Eventually, we will be interested in the limit when both $L$ and $\ell$ are sent to infinity, but $\ell/L \to 0$. We consider smooth field configurations $\psi$ with non-vanishing support in $\mathcal{I}_2$ and zero in the complement of $\mathcal{I}_1$: these field configurations can be seen alternatively as belonging either to the whole line problem by continuing the field as zero to infinity, or as the PBC problem by considering periodic repetitions of the fundamental interval. It is useful to conveniently switch between the two viewpoints. As long as the time evolution does not propagate the NLS field outside of $\mathcal{I}_1$, we can compare the whole-line transfer matrix (C3) with the periodic one, and in particular have $a(\lambda) = e^{iL\lambda/2} a^{\text{PBC}}(\lambda)$: boundary conditions on the PBC transfer matrix impose $a^{\text{PBC}}(\lambda_j) = \pm 1$, by taking the logarithm and using the representation (C5), one reaches [52]

$$\frac{I_j}{L} = \frac{\lambda_j}{2\pi} - \int \frac{d\lambda'}{2\pi} \frac{2}{\lambda_j - \lambda'} \frac{\log|a(\lambda')|}{L\pi}. \quad (C6)$$

The similarity with the logarithmic form of the Bethe equations (B3) is now evident, and this analogy, together with the fact that $\log|a(\lambda)|$ generates the conserved quantities, further encourages to identify the classical root density as $\rho(\lambda) \equiv \lim_{L \to \infty} \frac{1}{\pi} \log|a(\lambda)|$ [45, 52, 56], and the classical scattering phase $\Theta(\lambda) = \frac{2}{\lambda}$, then leading to the classical scattering kernel $\varphi(\lambda) = \partial_\lambda \Theta(\lambda)$. The argument presented here is not sufficient to derive the thermodynamics by paralleling the framework of TBA discussed in Appendix B, as Eq. (C6) is already in a coarse grained form and hinders the possibility of counting the microscopic solutions and characterize the classical entropy. A rigorous derivation passes through the identification of the action-angle microscopic variables [69], leading to the microscopic analogue of Eq. (C6). We will not go through this highly technical detour, as we provide an alternative derivation by considering the semiclassical limit of the LL in Appendix D. The final outcome is an entropy of the general form (B6), provided the classical root density is used and that the entropy factor is of Rayleigh-Jeans type $f(x) = \log x$.

## Appendix D: From Lieb-Liniger to Non Linear Schroedinger: semiclassical limits and PDF

For completeness, we provide a short summary of the semiclassical limit connecting the Lieb-Liniger model and the NLS, in the defocusing regime (see Ref. [54] for the focusing case). This Appendix summarizes the material presented in Ref. [45], to which we refer for a more in-depth discussion.

**The semiclassical limit of thermodynamics.**— The semiclassical limit is achieved for weak interactions and in the large occupation limit: in short, the energy scale should be such that the quantum model cannot resolve
the addition of a single excitation, which can be then coarse grained to a continuous variable. The semiclassical limit does not require integrability, and it is best introduced by looking at thermal expectation values. We consider a thermal partition function $\mathcal{Z} = \text{Tr}e^{-\beta_q(\hat{H}-\mu\hat{N})}$ with $\hat{N}$ the number of particles, and introduce a path integral representation of the quantum field $\hat{\psi}(x) \to \phi(\tau, x)$ with the Euclidean time having compact support $\tau \in [0, \beta_q]$, and conveniently rewrite the field in terms of the modes in the Euclidean time $\phi(\tau, x) = \sum_{n=-\infty}^{\infty} \frac{e^{i2\pi n\tau/\beta_q}}{\sqrt{c}} \psi_n(x)$. The semiclassical limit is achieved in a double limit of weak interaction $c \to 0$ and large temperature $\beta_q \to 0$, in such a way $\beta = \beta_q/c$, with $\beta$ the classical temperature. In terms of the partition function one has

$$ \mathcal{Z} = \int \mathcal{D}\phi e^{-\int dx \int_0^{\beta_q} d\tau \frac{1}{2}(\phi^* \partial_\tau \phi - \phi \partial_\tau \phi^*) + |\partial_x \phi|^2 + c|\phi|^4 + \mu|\phi|^2} = \int \prod_n \mathcal{D}\psi_n \exp \left\{ -\frac{\beta_q}{c} \int dx \sum_n \left[ \left( \frac{2\pi n}{\beta_q} \right)^2 + \mu \right] |\psi_n|^2 + |\partial_x \psi_n|^2 \right] + \sum_{n_1+n_2=n_3+n_4} \psi_{n_1}^* \psi_{n_2}^* \psi_{n_3} \psi_{n_4} \right\}. \quad (\text{D1}) $$

The double limit is immediately understood, as for small $\beta_q$ all the oscillating modes $\psi_{n \neq 0}$ are largely suppressed. The only non-vanishing field is the zero mode $\psi_{n=0}(x)$ which can be identified with the classical field of the NLS $\psi_{n=0}(x) \to \psi(x)$, and the quantum partition function becomes the classical partition function. This analysis can be extended to time evolution, and also to the integrability framework and GGEs with the correspondence $\beta_{q,n} \approx c\beta_n$. The mode decomposition over the Matsubara frequencies also establishes the correspondence between quantum and classical fields $\hat{\psi}(x) \approx \frac{1}{\sqrt{c}}\psi(x)$, which holds in a weak sense when computing correlation functions.

**The semiclassical limit of the TBA and PDF.**— The semiclassical limit can be taken at the level of thermodynamic Bethe ansatz, and the generating function for the moments of the density (B8) can be also projected to the classical limit finally leading to the formula for the PDF presented in the main text. We leave the detailed calculations to the Ref. [45, 52], and we limit ourselves to summarize the key-points. The correspondence at the leading order in the semiclassical limit between the quantum root density $\rho_q(\lambda)$, total root density $\rho_q^t(\lambda)$, and scattering kernel $\varphi_q(\lambda)$ with the classical counterparts is

$$ \begin{aligned} \rho_q(\lambda) &\approx \frac{1}{c}\rho(\lambda) \\ \rho_q^t(\lambda) - \rho_q(\lambda) &\approx \rho^t(\lambda) \\ \varphi_q(\lambda) &\approx \delta(\lambda) + c\varphi(\lambda) \end{aligned} \quad (\text{D2}) $$

All the general formulas of the quantum TBA becomes the analogue expression in the classical variables upon the above substitution. Furthermore, the semiclassical limit of the quantum entropy (B6) becomes the classical

entropy discussed at the end of Appendix C. The minimization of the free energy leads to the same integral equations (B7), provided the Rayleigh-Jeans distribution $f(x) = \log x$ is used. As discussed in Ref. [45], the semiclassical limit can be also extended to the formula for generating moments of $|\psi|^2$ (B8). However, this does not give directly the PDF: in Ref. [45] an equivalent representation of the moment generating function is used, which is easier to be inverted to finally obtain the sought formula for the PDF

$$ P_{\text{GGE}}(|\psi|^2 = d) = 16 \int_0^\infty d\theta J_0(8\sqrt{\theta d}) e^{-\frac{16}{\pi} \int_0^{\sqrt{\theta}} d\tau \tau S(\tau)}, \quad (\text{D3}) $$

with $J_0(x)$ a modified Bessel function of the first kind, the classical filling and total root density are $\vartheta(\lambda) = \rho(\lambda)/\rho^t(\lambda)$ and $\rho^t(\lambda) = \frac{1}{2\pi} - \int \frac{d\lambda'}{2\pi} \frac{2}{\lambda-\lambda'} \partial_{\lambda'}\rho(\lambda')$, and $S(\tau) \equiv \int d\lambda s_\tau(\lambda)$. The function $s_\tau(\lambda)$ is the solution of the auxiliary problem

$$ [\vartheta^{-1}(\lambda) + i0^{+}] s_{\tau}(\lambda) = 1 + \int \frac{d\lambda'}{2\pi} \frac{2}{\lambda - \lambda'} (4i\tau - \partial_{\lambda'}) s_{\tau}(\lambda'). \quad (\text{D4}) $$

The detailed derivation is left to the original reference [45]. Here, we provide a quick consistency check verifying that, in the limit of weak non-linearity, $P_{\text{GGE}}$. We achieve weak non-linearities in the limit of a broad plane-wave momentum distribution in the initial conditions $\Delta k \to \infty$. In this limit, the root density converges to the momentum distribution and we can therefore approximate it $\rho(\lambda) \approx \frac{1}{\Delta k \sqrt{\pi}} e^{-\lambda^2/\Delta k^2}$ and, at the leading order, $\rho^t(\lambda) \approx \frac{1}{2\pi}$, resulting in $\vartheta(\lambda) \approx \frac{2\sqrt{\pi}}{\Delta k} e^{-\lambda^2/\Delta k^2}$. We can now use this approximation in Eq. (D4) and notice that, since $\vartheta(\lambda)$ is very small, at the leading order in $\Delta k$ the solution for Eq. (D4) is $s_\tau(\lambda) \approx \vartheta(\lambda)$, leading to $S(\tau) \approx \frac{1}{2\pi}$. Importantly $S(\tau)$ is $\tau$-independent at this order. Using the approximate value of $S(\tau)$ in Eq. (D3) we obtain

$$ P_{\text{GGE}}(|\psi|^2 = d) \approx 16 \int_0^\infty d\theta J_0(8\sqrt{\theta d}) e^{-16\theta}, \quad (\text{D5}) $$

This last integral can be performed with Mathematica resulting in $P_{\text{GGE}}(|\psi|^2 = d) \approx e^{-d}$, as expected in the limit of zero non-linearity.

**How to compute the PDF in practice.**— We briefly discuss how to numerically solve the formula for the PDF: a working code is provided on Zenodo [66]. The TBA equations for computing the total root density and the auxiliary function $s_\tau(\lambda)$ can be discretized with standard methods and transformed in matrix-valued equations, see eg. [45] and [66]. The crucial point is correctly computing $e^{-\frac{16}{\pi} \int_0^{\sqrt{\theta}} d\tau S(\tau)}$: this function rapidly vanishes at large $\theta$ and thus allows for a truncation of the integral in Eq. (D3), but it develops zeroes associated with poles in $S(\tau)$. These poles need to be extracted and properly accounted for. We first properly rewrite (D4) in an operatorial form $\int d\lambda' [\Omega_\tau(\lambda, \lambda') + i0^+\delta(\lambda - \lambda')]s_\tau(\lambda') = 1$,
the definition of $\Omega_{\tau}(\lambda, \lambda')$ is readily obtained by a direct comparison with Eq. (D4). $\Omega_{\tau}(\lambda, \lambda')$ is parametrically dependent on $\tau$, and as an operator on the $\lambda$-space is Hermitian, and thus can be diagonalized. Let us call $b_{n;\tau}(\lambda)$ the eigenfunctions with eigenvalues $\mu_n(\tau)$, we assume $b_{n;\tau}(\lambda)$ are normalized to unity $\int d\lambda |b_{n;\tau}(\lambda)|^2$. In this representation, we can rewrite $S(\tau)$ as

$$S(\tau) = \sum_n \frac{B_n(\tau)}{\mu_n(\tau) + i0^+}, \quad B_n \equiv \int d\lambda b_{n;\tau}(\lambda). \quad (D6)$$

One needs now to isolate the possible singularities. Let us

fix a maximum cutoff $\tau < \sqrt{\theta_{\max}}$ for the computation of the PDF (D3). We tabulate the spectrum of the operator $\Omega_{\tau}(\lambda, \lambda')$ scanning through $\tau \in [0, \sqrt{\theta_{\max}}]$ and looking for solutions $\mu_n(\tau) = 0$, which are then tabulated as pairs $(n_j, \tau_j)$ in such a way $\mu_{n_j}(\tau_j) = 0$. We split $S(\tau)$ in a regular and singular part

$$\tau S(\tau) = \tau S^{\text{reg}}(\tau) + \sum_j \frac{\tau_j B_{n_j}(\tau_j)}{\partial_{\tau} \mu_{n_j}(\tau_j)(\tau - \tau_j) + i0^+}. \quad (D7)$$

As the next step, we can explicitly integrate over the poles when computing $\int d\tau \tau S(\tau)$. The explicit integration of one of these poles gives

$$\exp\left[-\frac{16}{\pi} \int_0^{\sqrt{\theta}} d\tau \, \tau \frac{\tau_j B_{n_j}(\tau_j)}{\partial_{\tau} \mu_{n_j}(\tau_j)(\tau - \tau_j) + i0^+}\right] = \left|\frac{\sqrt{\theta} - \tau_j}{\tau_j}\right|^{\frac{16}{\pi} \frac{\tau_j B_{n_j}(\tau_j)}{\partial_{\tau} \mu_{n_j}(\tau_j)}} \exp\left[\frac{16}{\pi} \frac{\tau_j B_{n_j}(\tau_j)}{|\partial_{\tau} \mu_{n_j}(\tau_j)|} i\pi H(\sqrt{\theta} - \tau_j)\right], \quad (D8)$$

with $H(x)$ the Heaviside Theta function $H(x > 0) = 1$ and zero otherwise. Apparently, the function above is non-analytic: however, it can be analytically shown [69] that $-\frac{16}{\pi} \frac{\tau_j B_{n_j}(\tau_j)}{\partial_{\tau} \mu_{n_j}(\tau_j)} = 1$, leading to simple zeroes $\exp\left[-\frac{16}{\pi} \int_0^{\sqrt{\theta}} d\tau \, \tau \frac{\tau_j B_{n_j}(\tau_j)}{\partial_{\tau} \mu_{n_j}(\tau_j)(\tau - \tau_j) + i0^+}\right] = \frac{\sqrt{\theta} - \tau_j}{\tau_j}$. Notice that in numerical discretizations of the kernel $\Omega_{\tau}(\lambda, \lambda')$, the identity $-\frac{16}{\pi} \frac{\tau_j B_{n_j}(\tau_j)}{\partial_{\tau} \mu_{n_j}(\tau_j)} = 1$ will not be exact, but will have small corrections due to the discretization. After having explicitly integrated over the singular terms, one

is left with

$$e^{-\frac{16}{\pi} \int_0^{\sqrt{\theta}} d\tau \, \tau S(\tau)} = e^{-\frac{16}{\pi} \int_0^{\sqrt{\theta}} d\tau \, \tau S^{\text{reg}}(\tau)} \prod_j \left( \frac{\sqrt{\theta} - \tau_j}{\tau_j} \right), \quad (D9)$$

This is the most convenient expression for a numerical implementation, where in the first step one finds the zeroes in the spectrum of $\Omega_{\tau}$ finding $\tau_j$ and removing the singularities from $S(\tau)$, then the integral over the regular part of $S^{\text{reg}}(\tau)$ can be safely evaluated. This implementation is used in the code provided on Zenodo [66].

[1] E. T. Jaynes, Information theory and statistical mechanics, Phys. Rev. 106, 620 (1957).

[2] E. T. Jaynes, Information theory and statistical mechanics. ii, Phys. Rev. 108, 171 (1957).

[3] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).

[4] M. Serbyn, D. A. Abanin, and Z. Papić, Quantum many-body scars and weak breaking of ergodicity, Nature Physics 17, 675 (2021).

[5] A. Gromov and L. Radzihovsky, Colloquium: Fracton matter, Rev. Mod. Phys. 96, 011001 (2024).

[6] I. Bloch, J. Dalibard, and S. Nascimbène, Quantum simulations with ultracold quantum gases, Nature Physics 8, 267 (2012).

[7] X.-W. Guan and P. He, New trends in quantum integrability: recent experiments with ultracold atoms, Reports on Progress in Physics 85, 114001 (2022).

[8] P. Suret, S. Randoux, A. Gelash, D. Agafontsev, B. Doyon, and G. El, Soliton gas: Theory, numerics, and experiments, Phys. Rev. E 109, 061001 (2024).

[9] F. A. Smirnov, Form Factors in Completely Integrable Models of Quantum Field

Theory (WORLD SCIENTIFIC, 1992) https://www.worldscientific.com/doi/pdf/10.1142/1115.

[10] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, Quantum Inverse Scattering Method and Correlation Functions, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 1993).

[11] L. D. Faddeev and L. A. Takhtajan, Hamitonian methods in the theory of solitons (1987).

[12] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).

[13] S. Sotiriadis and P. Calabrese, Validity of the gge for quantum quenches from interacting to noninteracting models, Journal of Statistical Mechanics: Theory and Experiment 2014, P07024 (2014).

[14] B. Doyon, Thermalization and pseudolocality in extended quantum systems, Communications in Mathematical Physics 351, 155 (2017).

[15] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a completely integrable many-body quantum system: An ab initio study of the dynamics of the highly excited states of 1d lattice hard-core bosons, Phys. Rev.
Lett. 98, 050405 (2007).

[16] P. Calabrese, F. H. L. Essler, and M. Fagotti, Quantum quench in the transverse-field ising chain, Phys. Rev. Lett. 106, 227203 (2011).

[17] M. Kormos, M. Collura, and P. Calabrese, Analytic results for a quantum quench from free to hard-core one-dimensional bosons, Phys. Rev. A 89, 013609 (2014).

[18] B. Pozsgay, M. Mestyán, M. A. Werner, M. Kormos, G. Zaránd, and G. Takács, Correlations after quantum quenches in the xxz spin chain: Failure of the generalized gibbs ensemble, Phys. Rev. Lett. 113, 117203 (2014).

[19] F. H. L. Essler and M. Fagotti, Quench dynamics and relaxation in isolated integrable quantum spin chains, Journal of Statistical Mechanics: Theory and Experiment 2016, 064002 (2016).

[20] E. Ilievski, M. Medenjak, T. Prosen, and L. Zadnik, Quasilocal charges in integrable lattice systems, Journal of Statistical Mechanics: Theory and Experiment 2016, 064008 (2016).

[21] P. Calabrese, F. H. L. Essler, and G. Mussardo, Introduction to ‘quantum integrability in out of equilibrium systems’, Journal of Statistical Mechanics: Theory and Experiment 2016, 064001 (2016).

[22] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer, Experimental observation of a generalized gibbs ensemble, Science 348, 207 (2015), https://www.science.org/doi/pdf/10.1126/science.1257026.

[23] M. Takahashi, Thermodynamics of one-dimensional solvable models (Cambridge University Press, 2005).

[24] E. Ilievski, E. Quinn, J. D. Nardis, and M. Brockmann, String-charge duality in integrable lattice models, Journal of Statistical Mechanics: Theory and Experiment 2016, 063101 (2016).

[25] J. M. Wilson, N. Malvania, Y. Le, Y. Zhang, M. Rigol, and D. S. Weiss, Observation of dynamical fermionization, Science 367, 1461 (2020), https://www.science.org/doi/pdf/10.1126/science.aaz0242.

[26] N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, Generalized hydrodynamics in strongly interacting 1d bose gases, Science 373, 1129 (2021), https://www.science.org/doi/pdf/10.1126/science.abf0147.

[27] L. Dubois, G. Thémèze, F. Nogrette, J. Dubail, and I. Bouchoule, Probing the local rapidity distribution of a one-dimensional bose gas, Phys. Rev. Lett. 133, 113402 (2024).

[28] Y. Le, Y. Zhang, S. Gopalakrishnan, M. Rigol, and D. S. Weiss, Observation of hydrodynamization and local prethermalization in 1d bose gases, Nature 618, 494 (2023).

[29] K.-Y. Li, Y. Zhang, K. Yang, K.-Y. Lin, S. Gopalakrishnan, M. Rigol, and B. L. Lev, Rapidity and momentum distributions of one-dimensional dipolar quantum gases, Phys. Rev. A 107, L061302 (2023).

[30] M. Horvath, A. Bastianello, S. Dhar, R. Koch, Y. Guo, J.-S. Caux, M. Landini, and H.-C. Nägerl, Observing bethe strings in an attractive bose gas far from equilibrium (2025), arXiv:2505.10550 [cond-mat.quant-gas].

[31] V. E. Zakharov, Turbulence in integrable systems, Studies in Applied Mathematics 122, 219 (2009), https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1467-9590.2009.00430.x.

[32] D. Agafontsev and V. E. Zakharov, Integrable turbulence and formation of rogue waves, Nonlinearity 28, 2791 (2015).

[33] A. Gelash, D. Agafontsev, V. Zakharov, G. El, S. Randoux, and P. Suret, Bound state soliton gas dynamics underlying the spontaneous modulational instability, Physical review letters 123, 234102 (2019).

[34] T. Congy, G. A. El, G. Roberti, A. Tovbis, S. Randoux, and P. Suret, Statistics of extreme events in integrable turbulence, Phys. Rev. Lett. 132, 207201 (2024).

[35] A. Tikan, S. Bielawski, C. Szwaj, S. Randoux, and P. Suret, Single-shot measurement of phase and amplitude by using a heterodyne time-lens system and ultrafast digital time-holography, Nature Photonics 12, 228 (2018).

[36] A. E. Kraych, D. Agafontsev, S. Randoux, and P. Suret, Statistical properties of the nonlinear stage of modulation instability in fiber optics, Phys. Rev. Lett. 123, 093902 (2019).

[37] G. Michel, F. Bonnefoy, G. Ducrozet, G. Prabhudesai, A. Cazaubiel, F. Copie, A. Tikan, P. Suret, S. Randoux, and E. Falcon, Emergence of peregrine solitons in integrable turbulence of deep water gravity waves, Phys. Rev. Fluids 5, 082801 (2020).

[38] I. Redor, H. Michallet, N. Mordant, and E. Barthélemy, Experimental study of integrable turbulence in shallow water, Phys. Rev. Fluids 6, 124801 (2021).

[39] P. Suret, A. Tikan, F. Bonnefoy, F. m. c. Copie, G. Ducrozet, A. Gelash, G. Prabhudesai, G. Michel, A. Cazaubiel, E. Falcon, G. El, and S. Randoux, Nonlinear spectral synthesis of soliton gas in deep-water surface gravity waves, Phys. Rev. Lett. 125, 264101 (2020).

[40] S. Randoux, P. Walczak, M. Onorato, and P. Suret, Intermittency in integrable turbulence, Phys. Rev. Lett. 113, 113902 (2014).

[41] S. Randoux, P. Walczak, M. Onorato, and P. Suret, Nonlinear random optical waves: Integrable turbulence, rogue waves and intermittency, Physica D 333, 323 (2016).

[42] G. P. Agrawal, Nonlinear Fiber Optics, fifth edition ed. (Elsevier/Academic Press, Amsterdam, 2013).

[43] Y. Bromberg, Y. Lahini, E. Small, and Y. Silberberg, Hanbury brown and twiss interferometry with interacting photons, Nature Photonics 4, 721 (2010).

[44] S. Derevyanko and E. Small, Nonlinear propagation of an optical speckle field, Phys. Rev. A 85, 053816 (2012).

[45] G. D. V. D. Vecchio, A. Bastianello, A. D. Luca, and G. Mussardo, Exact out-of-equilibrium steady states in the semiclassical limit of the interacting Bose gas, SciPost Phys. 9, 002 (2020).

[46] T. Bonnemain, B. Doyon, and G. El, Generalized hydrodynamics of the kdv soliton gas, Journal of Physics A: Mathematical and Theoretical 55, 374004 (2022).

[47] V. Zakharov, Kinetic equation for solitons, Sov. Phys. JETP 33, 538 (1971).

[48] G. A. El and A. M. Kamchatnov, Kinetic equation for a dense soliton gas, Phys. Rev. Lett. 95, 204101 (2005).

[49] G. A. El, Soliton gas in integrable dispersive hydrodynamics, Journal of Statistical Mechanics: Theory and Experiment 2021, 114001 (2021).

[50] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Emergent hydrodynamics in integrable quantum systems out of equilibrium, Phys. Rev. X 6, 041065 (2016).

[51] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Transport in out-of-equilibrium xxz chains: Exact profiles of charges and currents, Phys. Rev. Lett. 117, 2791 (2016).
207201 (2016).

[52] A. D. Luca and G. Mussardo, Equilibration properties of classical integrable field theories, Journal of Statistical Mechanics: Theory and Experiment **2016**, 064011 (2016).

[53] A. Bastianello, B. Doyon, G. Watts, and T. Yoshimura, Generalized hydrodynamics of classical integrable field theory: the sinh-Gordon model, SciPost Phys. **4**, 045 (2018).

[54] R. Koch, J.-S. Caux, and A. Bastianello, Generalized hydrodynamics of the attractive non-linear schrödinger equation, Journal of Physics A: Mathematical and Theoretical **55**, 134001 (2022).

[55] R. Koch and A. Bastianello, Exact thermodynamics and transport in the classical sine-Gordon model, SciPost Phys. **15**, 140 (2023).

[56] Y. Bezzaz, L. Dubois, and I. Bouchoule, Rapidity distribution within the defocusing non-linear Schrödinger equation model, SciPost Phys. Core **6**, 064 (2023).

[57] A. Bastianello, Z. Krajnik, and E. Ilievski, Landau-lifschitz magnets: Exact thermodynamics and transport, Phys. Rev. Lett. **133**, 107102 (2024).

[58] P. Suret, R. E. Koussaifi, A. Tikan, C. Evain, S. Randoux, C. Szwaj, and S. Bielawski, Single-shot observation of optical rogue waves in integrable turbulence using time microscopy, Nature Communications **7**, 13136 (2016).

[59] A. Lebel, A. Tikan, S. Randoux, P. Suret, and F. Copie, Single-shot observation of breathers from noise-induced modulation instability using heterodyne temporal imaging, Optics Letters **46**, 298 (2021).

[60] P. Dorey, Exact s-matrices, in *Conformal Field Theories and Integrable Models*, edited by Z. Horváth and L. Palla (Springer Berlin Heidelberg, Berlin, Heidelberg, 1997) pp. 85-125.

[61] C. Boldrighini, R. L. Dobrushin, and Y. M. Sukhov, One-dimensional hard rod caricature of hydrodynamics, Journal of Statistical Physics **31**, 577 (1983).

[62] B. Doyon and H. Spohn, Dynamics of hard rods with initial domain wall state, Journal of Statistical Mechanics: Theory and Experiment **2017**, 073210 (2017).

[63] B. Doyon, T. Yoshimura, and J.-S. Caux, Soliton gases and generalized hydrodynamics, Phys. Rev. Lett. **120**, 045301 (2018).

[64] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler, and T. Prosen, Complete generalized gibbs ensembles in an interacting theory, Phys. Rev. Lett. **115**, 157201 (2015).

[65] J. Yang, *Nonlinear waves in integrable and nonintegrable systems* (SIAM, 2010).

[66] A. Bastianello, A. Tikan, F. Copie, S. Randoux, and P. Suret, Observation of a generalized Gibbs ensemble in photonics, Zenodo (2025).

[67] G. A. El and M. A. Hoefer, Dispersive shock waves and modulation theory, Physica D: Nonlinear Phenomena Dispersive Hydrodynamics, **333**, 11 (2016).

[68] Y. S. Kivshar and B. Luther-Davies, Dark optical solitons: Physics and applications, Physics Reports **298**, 81 (1998).

[69] R. K. Bullough, D. J. Pilling, and J. Timonen, Quantum and classical statistical mechanics of the sinh-gordon equation, Journal of Physics A: Mathematical and General **19**, L955 (1986).

[70] A. Tovbis and F. Wang, Spectral theory of soliton gases for the defocusing nls equation (2025), arXiv:2503.01132 [nlin.PS].

[71] A. Bastianello, L. Piroli, and P. Calabrese, Exact local correlations and full counting statistics for arbitrary states of the one-dimensional interacting bose gas, Phys. Rev. Lett. **120**, 190601 (2018).

[72] A. Bastianello and L. Piroli, From the sinh-gordon field theory to the one-dimensional bose gas: exact local correlations and full counting statistics, Journal of Statistical Mechanics: Theory and Experiment **2018**, 113104 (2018).

[73] E. H. Lieb and W. Liniger, Exact analysis of an interacting bose gas. i. the general solution and the ground state, Phys. Rev. **130**, 1605 (1963).

[74] E. H. Lieb, Exact analysis of an interacting bose gas. ii. the excitation spectrum, Phys. Rev. **130**, 1616 (1963).

[75] T. Kinoshita, T. Wenger, and D. S. Weiss, A quantum newton's cradle, Nature **440**, 900 (2006).

[76] J. Esteve, J.-B. Trebbia, T. Schumm, A. Aspect, C. I. Westbrook, and I. Bouchoule, Observations of density fluctuations in an elongated bose gas: Ideal gas and quasicondensate regimes, Phys. Rev. Lett. **96**, 130403 (2006).

[77] J. Armijo, T. Jacqmin, K. V. Kheruntsyan, and I. Bouchoule, Probing three-body correlations in a quantum gas using the measurement of the third moment of density fluctuations, Phys. Rev. Lett. **105**, 230402 (2010).

[78] B. Doyon, Exact large-scale correlations in integrable systems out of equilibrium, SciPost Phys. **5**, 054 (2018).

[79] J. De Nardis, B. Doyon, M. Medenjak, and M. Panfil, Correlation functions and transport coefficients in generalized hydrodynamics, Journal of Statistical Mechanics: Theory and Experiment **2022**, 014002 (2022).

[80] F. Franchini, *An introduction to integrable techniques for one-dimensional quantum systems*, Vol. 940 (Springer, 2017) arXiv:1609.02100 [cond-mat.stat-mech].

[81] H. Bethe, Zur theorie der metalle, Zeitschrift für Physik **71**, 205 (1931).

[82] M. G. Forest and D. W. McLaughlin, Spectral theory for the periodic sine-gordon equation: A concrete viewpoint, Journal of Mathematical Physics **23**, 1248 (1982), https://pubs.aip.org/aip/jmp/article-pdf/23/7/1248/19158407/1248_1_online.pdf.