# Spatio-temporal thermalization and adiabatic cooling of guided light waves

Lucas Zanaglia, Josselin Garnier, Iacopo Carusotto, Valérie Doya, Claire Michel, Antonio Picozzi

► To cite this version:

Lucas Zanaglia, Josselin Garnier, Iacopo Carusotto, Valérie Doya, Claire Michel, et al. Spatio-temporal thermalization and adiabatic cooling of guided light waves. 2025. hal-05138789

HAL Id: hal-05138789

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

Preprint submitted on 1 Jul 2025

**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 Attribution 4.0 International License
Spatio-temporal thermalization and adiabatic cooling of guided light waves

Lucas Zanaglia¹, Josselin Garnier², Iacopo Carusotto³, Valérie Doya¹, Claire Michel¹,⁴, Antonio Picozzi⁵

¹ Université Côte d'Azur, CNRS, Institut de Physique de Nice, Nice, France

² CMAP, CNRS, Ecole polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

³Pitaevskii BEC Center, CNR-INO and Dipartimento di Fisica, Università di Trento, I-38123 Trento, Italy.

⁴Institut Universitaire de France (IUF), 1 rue Descartes, 75005 Paris, France and

⁵ Université Bourgogne Europe, CNRS, Laboratoire Interdisciplinaire Carnot de Bourgogne ICB UMR 6303, 21000 Dijon, France

We propose and theoretically characterize three-dimensional spatio-temporal thermalization of a continuous-wave classical light beam propagating along a multi-mode optical waveguide. By combining a non-equilibrium kinetic approach based on the wave turbulence theory and numerical simulations of the field equations, we anticipate that thermalizing scattering events are dramatically accelerated by the combination of strong transverse confinement with the continuous nature of the temporal degrees of freedom. In connection with the blackbody catastrophe, the thermalization of the classical field in the continuous temporal direction provides an intrinsic mechanism for adiabatic cooling and, then, spatial beam condensation. Our results open new avenues in the direction of a simultaneous spatial and temporal beam cleaning.

*Introduction.*- Fluids of light in propagating geometries are an emerging platform to study the physics of quantum gases in novel regimes [1, 2]. A hydrodynamic reformulation of the two-dimensional (2D) transverse dynamics of monochromatic light beams underlies the observation of a variety of physical effects, such as complex vortex dynamics [3], superfluid light behaviours [4, 5], the formation of prethermalized equilibrium states [6], and even optical analogs of gravitational phenomena [7–9]. On the basis of spatial beam cleaning experiments [10–12] and related theoretical works [13–19], the relaxation to a fully thermalized equilibrium state has recently been experimentally observed for monochromatic light propagating along parabolic-shaped (graded-index) multimode optical fibers [20–23]. This intense activity has spurred the emerging key area of optical thermodynamics [13], marked by notable advances, such as negative temperature states [13, 24–26], calorimetry of photon gases [27], and their Joule-Thomson expansion [28], or the development of nonequilibrium approaches [29–32].

The physics is far richer when one goes beyond the monochromatic assumption and allows for a spatio-temporal 3D dynamics. In this case, the mapping of light propagation onto a fluid of light can be rigorously formulated at the microscopic quantum level leading to a full quantum fluid theory of light [33]. In synergy with pioneering experiments on quantum fluctuation features [34], the study of hydrodynamic effects in 3D fluids of light is closely intertwined with the strong on-going activity on the spatio-temporal (ST) dynamics of light beams, either in bulk media (X-waves [35], ST solitons [36], ST self-similarity [37]), or in multimode waveguides [38] (multimode solitons [39–41], conical emission [42–44], instabilities [45, 46], supercontinuum generation [10, 47, 48]).

One of the most exciting challenges in this context is the observation of full ST thermalization of light, with the consequent intriguing possibility of a spectral condensation also in the time domain, in combination with

the transverse spatial beam-cleaning. While ST thermalization of quantum light is expected to require impractically long propagation distances [49], general theoretical results on the 3D thermalization of classical waves [50–54] suggest that thermalization may be strongly accelerated in the case of classical light waves by bosonic stimulation effects. A most serious obstacle in directly transposing these results to the optical case is spatial diffraction which, in a spatially homogeneous (bulk) medium, leads to a quick expansion of any realistic light beam and effectively blocks the thermalization process [55].

In this Letter, taking inspiration from on-going efforts in the cold-atom context [56–60], we propose a strategy to overcome these difficulties by considering a multimode waveguide geometry, in which we theoretically anticipate an efficient 3D thermalization of classical waves in both the spatial and the temporal directions. In such systems, the transverse dynamics is confined in space by a suitable refractive index profile, which prevents the beam from expanding and, thanks to a strong spatial overlap of the modes, guarantees efficient nonlinear interactions. On the other hand, by working with a temporally continuous incoherent input wave, rather than the usual coherent pulsed light [16, 20–23, 25–29], the light beam remains homogeneous in the temporal direction, with no need for confinement along this axis to overcome pulse broadening effects. As a result, while the considered transverse confinement does not provide resonances allowing for a solely spatial thermalization, ST scattering events have a continuum of temporal modes at their disposal, and thus do not suffer from the obstacles found in the thermalization of systems with a discrete spectrum [53, 61, 62]. This educated guess is confirmed by our calculations: While the purely spatial system does not exhibit thermalization, the corresponding ST system exhibits a fast convergence toward a local Rayleigh-Jeans ST equilibrium, which is shown to drive an adiabatic cooling and spatial beam condensation of the incoherent field.
**The UPE model.** We present the theory of ST thermalization in the framework of a general model of light propagation that goes beyond the usual slowly-varying envelope and paraxial approximations inherent to nonlinear Schrödinger-like equations (NLSE). The model, known as unidirectional propagation equation (UPE), has been introduced with the aim of describing extreme nonlinear optics, down to the optical cycle time scale and the spatial scale of the wavelength of light [63–67]. The starting point are the Maxwell equations with the weak assumptions of scalar approximation and unidirectional forward light propagation (along $z > 0$). By expanding the electric field on the eigenmodes of the waveguide at central frequency $\omega_o$, we derive a modified version of the UPE governing the $z$-evolution of the modal amplitudes $b_m(\omega, z)$ (see Appendix I):

$$i\partial_z b_m = (\tilde{\beta}_m(\omega) + \omega/v_g)b_m - \gamma\Gamma_m(\omega)P_m(\omega, \mathbf{b}). \quad (1)$$

The propagation constants $\tilde{\beta}_m(\omega)$ describing linear dispersion effects, as well as the nonlinear dispersion coefficients $\Gamma_m(\omega)$, are functions of the frequency $\omega$:

$$\tilde{\beta}_m(\omega) = k_o - \sqrt{k^2(\omega) - 2k_o\beta_m}, \quad (2)$$

$$\Gamma_m(\omega) = \frac{\omega^2 k_o}{\omega_o^2 \sqrt{k^2(\omega) - 2k_o\beta_m}}, \quad (3)$$

where $k(\omega) = n(\omega)\omega/c$ with $n(\omega)$ being the refractive index and $c$ the vacuum speed of light, $k_o = k(\omega_o)$, $\beta_m = \tilde{\beta}_m(\omega_o)$, and $\gamma$ the nonlinear coefficient. In Fig. 1(b) (inset) we report a typical frequency dependence of $\tilde{\beta}_m(\omega)$ (with a central frequency $\omega_o$ located in the anomalous dispersion regime). The UPE is written in a reference frame that propagates with inverse group-velocity $v_g^{-1} = \partial_\omega k(\omega_o)$. We consider a cubic (Kerr) nonlinearity $P_m(\omega, \mathbf{b}) = \sum_{pqr} W_{mpqr} \int b_p(\omega_1, z)b_q^*(\omega_2, z)b_r(\omega_3, z)\delta_\omega d\omega_1 d\omega_2 d\omega_3$, where $W_{mpqr}$ accounts for the spatial overlap among the waveguide eigenmodes and $\delta_\omega = \delta(\omega - \omega_1 + \omega_2 - \omega_3)$.

The UPE is the most general scalar propagation equation that encompasses different approximate models of light propagation [65]. By expanding $\tilde{\beta}_m(\omega)$ and $\Gamma_m(\omega)$ around $\omega_o$, one can derive from the UPE the forward Maxwell equation, the first-order propagation equation, the nonlinear envelope equation (NEE), or the NLSE [65].

**Resonances induced by ST coupling.** Nonequilibrium light thermalization crucially relies on efficient quasi-resonances among quartets of modes [15]. For instance, observation of spatial thermalization has been reported in parabolic-shaped multimode fibers where such resonances are naturally present (see the review [12]), and can be suppressed in multimode step-index waveguides [15]. In what follows we show that in the ST case efficient quasi-resonances are recovered in generic waveguide configurations, leading to efficient ST thermalization. By considering a step-index waveguide, we report in Fig. 1(a) the histogram of four-mode resonances for

<figure><img src="image_2.png" /><figcaption>FIG. 1: **Quasi-resonances.** Schematic visualization and corresponding histograms of four-mode resonances in (a) the pure (monochromatic) spatial case, and in (b) the (broadband polychromatic) ST case. Owing to the additional temporal degrees of freedom, in the ST case the histogram collapses to efficient quasi-resonances, $|\Delta\beta^{ST}L_{nl}| \ll 1$ (note the different vertical scales). Inset: Modal dispersion relations $\tilde{\beta}_m(\omega)$ (in $\text{mm}^{-1}$) from Eq.(2) for a step-index waveguide supporting 26 modes.</figcaption></figure>

the (monochromatic) spatial case, $\Delta\beta^S = \Delta\beta_{mpqr} = \beta_m + \beta_q - \beta_p - \beta_r$. Clearly, the occurrence of four-modes quasi-resonances verifying $|\Delta\beta_{mpqr}| \ll 1/L_{nl}$ is deficient, where $L_{nl} \simeq 1/(|\gamma|\tilde{N})$ is the nonlinear length, and $\tilde{N}$ the average power (see Appendix I). The introduction of the temporal degrees of freedom fundamentally changes the physical picture, because the propagation constants become functions of the frequency, $\beta_m \to \tilde{\beta}_m(\omega)$, see the inset of Fig. 1(b). In this broadband polychromatic case, the system has an extra degree of freedom to fulfill efficient resonances $|\Delta\beta^{ST}L_{nl}| \ll 1$, where $\Delta\beta^{ST} = \min(\Delta\tilde{\beta}_{mpqr}^{123})$ is the most favorable resonance (with the minimum taken over the frequencies $\{\omega, \omega_{1,2,3}\}$), and $\Delta\tilde{\beta}_{mpqr}^{123} = \tilde{\beta}_m(\omega) - \tilde{\beta}_p(\omega_1) + \tilde{\beta}_q(\omega_2) - \tilde{\beta}_r(\omega_3)$. In contrast with the pure spatial case, the ST histogram collapses to nearly exact resonances, as remarkably shown in Fig. 1(b). The continuous time variable then enables the system to find exact resonances that offset the large frequency mismatch $\Delta\beta^S$ of purely spatial resonances.

**Spatio-temporal wave turbulence.** The intuitive physical picture about resonant interactions depicted in Fig. 1 can be formalized in the framework of the wave turbulence theory, which provides a detailed nonequilibrium description of the irreversible thermalization process [53, 62, 68–75]. In the purely spatial case, as anticipated in Fig. 1,
<figure><img src="image_3.png" /><figcaption>**FIG. 2:** Spatio-temporal thermalization. (a) Simulation of Eq.(1) in the NEE approximation: Evolution of the spatial spectrum $N_m^{ST}(z)$, showing the relaxation to the equilibrium RJ distribution Eq.(6) (dashed black line). (b) Evolution of the distance $D(z)$ to equilibrium, whose decrease to zero evidences ST thermalization. (c) This irreversible process is also characterized by a monotonous growth of entropy $S(z)$, as described by the H-theorem of the wave turbulence kinetic Eq.(4). The distance and entropy evolutions are in contrast with those of the spatial case, see Figs. 3(b)-(c). Temporal spectrum $|b_m|^2(\omega, z)$ of the fundamental mode $m = 0$ (d), intermediate mode $m = 10$ (e), highest mode ($m = 25$) (f), at $z = 0$ (dark blue) and $z = 570L_{nl}$ (red), showing thermalization to the RJ spectra (dashed light blue). Parameters: step-index waveguide supporting 26 modes (see the inset of Fig. 1(b)), with anomalous dispersion and defocusing nonlinearity, $\tau_0 = \sqrt{|\kappa_2|L_{nl}/2}$, $L_{nl} = 0.4m$, $\tilde{\omega}_c = 40/\tau_0$.</figcaption></figure>

the discrete set of resonant interactions suppresses quasi-resonances and ultimately freezes the thermalization process. This aspect was discussed in the context of discrete wave turbulence [53, 61], e.g., through the analysis of Fermi-Pasta-Ulam chains [62], or nonlinear disordered systems [76–80]. Accordingly, the discrete nature of wave turbulence can prevent the derivation of a continuous kinetic equation describing spatial-only thermalization, see Appendix III. Here, we show that the continuous nature of the temporal degrees of freedom restores efficient ST resonances, enabling the derivation of a hybrid discrete-continuous ST kinetic equation, involving discrete sums over the spatial modes and continuous integrals over the temporal spectrum.

It proves convenient to symmetrize the UPE by introducing the renormalized modal amplitudes $\tilde{b}_m(\omega) = b_m(\omega)/\sqrt{\Gamma_m(\omega)}$. In the weakly nonlinear regime relevant to usual multimode optical fiber experiments [12, 20–

<figure><img src="image_4.png" /><figcaption>**FIG. 3:** Pure spatial dynamics: Frozen thermalization. (a) Simulation of Eq.(29): Evolution of the spatial spectrum $N_m^S(z)$ starting from the same initial condition as in the ST simulation in Fig. 2. The thermalization process is frozen, as evidenced by the distance $D(z)$ to RJ equilibrium (b), and the entropy (c), whose evolutions are in contrast to the ST case in Fig. 2(b)-(c). Because of the large fluctuations of individual realizations, an average has been taken over 21 realizations.</figcaption></figure>

22, 25, 27], we derive the kinetic equation for the evolution of the ST spectrum $\langle \tilde{b}_m(\omega) \tilde{b}_p^*(\omega') \rangle = \tilde{n}_m(\omega) \delta_{mp} \delta(\omega - \omega')$ [$\delta_{mp}$ being the Kronecker symbol] (see Appendix II):

$$ \partial_z \tilde{n}_m(\omega) = 4\pi\gamma^2 \sum_{pqr} \int d\omega_1-3 |L_{mpqr}^{\omega_{123}}|^2 \mathbf{M}_{mpqr}(\tilde{n}) \times \delta(\omega - \omega_1 + \omega_2 - \omega_3) \delta(\Delta\tilde{\beta}_{mpqr}^{\omega_{123}}) \quad (4) $$

with the nonlinear interaction tensor

$$ L_{mpqr}^{\omega_{123}} = W_{mpqr} \sqrt{\Gamma_m(\omega) \Gamma_p(\omega_1) \Gamma_q(\omega_2) \Gamma_r(\omega_3)} \quad (5) $$

and the cubic nonlinear term $M_{mpqr}(\tilde{n}) = \tilde{n}_p(\omega_1)\tilde{n}_q(\omega_2)\tilde{n}_r(\omega_3) + \tilde{n}_m(\omega)\tilde{n}_p(\omega_1)\tilde{n}_r(\omega_3) - \tilde{n}_m(\omega)\tilde{n}_q(\omega_2)\tilde{n}_r(\omega_3) - \tilde{n}_m(\omega)\tilde{n}_p(\omega_1)\tilde{n}_q(\omega_2)$. Eq.(4) conserves the power $N_{\tilde{n}} = \sum_m \int \tilde{n}_m(\omega) d\omega$, the momentum $P_{\tilde{n}} = \sum_m \int \omega \tilde{n}_m(\omega) d\omega$, and the kinetic energy $E_{\tilde{n}} = \sum_m \int \tilde{\beta}_m(\omega) \tilde{n}_m(\omega) d\omega$. It exhibits a H-theorem of entropy growth $\partial_z S_{\tilde{n}}(z) \ge 0$, for the nonequilibrium entropy $S_{\tilde{n}}(z) = \sum_m \int \log(\tilde{n}_m(\omega)) d\omega$. Then at variance with the UPE (1) that is formally reversible (Hamiltonian) in the 'time' variable $z$, the kinetic Eq.(4) describes the actual nonequilibrium process of ST thermalization toward equilibrium. In terms of the original variables $b_m$, the equilibrium ST distribution has a Rayleigh-Jeans (RJ) form

$$ n_m^{\text{RJ}}(\omega) = \tilde{n}_m^{\text{RJ}}(\omega) \Gamma_m(\omega) = \frac{T \Gamma_m(\omega)}{\tilde{\beta}_m(\omega) + \omega(1/v_g - \lambda) - \mu}, \quad (6) $$

where the temperature $T$, chemical potential $\mu$, and average 'velocity' $\lambda$, are determined from the three conserved quantities ($N_{\tilde{n}}$, $P_{\tilde{n}}$, $E_{\tilde{n}}$).
<figure><img src="image_5.png" /><figcaption>FIG. 4: Local-equilibrium route to ST thermalization and adiabatic cooling. (a) Mode-integrated temporal spectrum of the field $\sum_m |b_m(\omega)|^2$ at $z = 160L_{nl}$ (blue), and local RJ equilibrium distribution over the reduced frequency window $[\tilde{\omega}_{c,-}^{loc}, \tilde{\omega}_{c,+}^{loc}]$ (orange). (b) Modal population $N_m^{loc}(z)/N$ computed from the local RJ equilibrium (dashed lines), and modal population $N_m(z)/N$ in the NEE simulation of Fig. 2 (continuous lines). (c) Evolution during propagation of the positive $\tilde{\omega}_{c,+}^{loc}(z)$, and negative $\tilde{\omega}_{c,-}^{loc}(z)$, frequency cut-off. (d) Local temperature $T^{loc}(z)$ as obtained from a fit of the numerical distribution with a RJ distribution within the frequency window $[\tilde{\omega}_{c,-}^{loc}, \tilde{\omega}_{c,+}^{loc}]$: The decrease in $T^{loc}(z)$ reflects an adiabatic cooling, which drives a spatial beam condensation characterized by the growth of $N_0^{loc}(z)/N$ in (c).</figcaption></figure>

**Spatio-temporal simulations.** To put these ideas on solid grounds, we have performed ST simulations by considering the example of the NEE model, because it provides the minimal ingredients for a non-trivial space-time coupling, whereby different frequency components (colours) diffract along different cone angles. The NEE is obtained from the UPE (1) with the approximations: $\tilde{\beta}_m(\tilde{\omega}) = \beta_m(1+\tilde{\omega}/\omega_0)^{-1} - \sum_{j \ge 2} \kappa_j \tilde{\omega}^2$, where $\kappa_j = \frac{1}{j!} \partial_j k(\omega_0)$ denote high-order dispersion effects, and $\Gamma_m(\tilde{\omega}) = 1+\tilde{\omega}/\omega_0$, with $\tilde{\omega} = \omega - \omega_0$ the frequency offset (see Appendix I). At variance with the experiments carried out so far to study light thermalization [12, 16, 20–23, 25–29], here the injected beam is spatio-temporally incoherent: the different $(\omega, m)$-components $b_m(\omega, z=0)$ are independent complex-valued Gaussian random variables of zero mean; each spatial mode $m$ has a Gaussian spectrum as a function of $\omega$ with the same FWHM and different amplitudes, see dark blue curves in Fig. 2(a,d-f).

Spatio-temporal thermalization to the Rayleigh-Jeans equilibrium (6) is illustrated in Fig. 2. Here, the parameters $(T, \lambda, \mu)$ in Eq.(6) are computed by considering the frequency cutoff $\tilde{\omega}_c \tau_0 = 40$ of the spectral grid used in the simulation (see Appendix IV). We report in Fig. 2(a) the evolution of the spatial mode distribution by integrating over the temporal frequencies, $N_m^{ST}(z) = \frac{1}{2\pi} \int |b_m(\omega, z)|^2 d\omega$. ST thermalization is confirmed

by the evolution of the distance to RJ equilibrium, $D(z) = \sum_p |N_p(z) - N_p^{RJ}|/\sum_p (N_p(z) + N_p^{RJ})$, which decreases to zero during propagation (note that $D$ is bounded, $0 \le D \le 1$), while the temporal spectra tend to converge to those predicted by the RJ equilibrium, see Fig. 2(d-f). Note that, to avoid the formation of temporal solitons that can freeze the thermalization process, in the anomalous dispersion regime of Fig. 2 a defocusing nonlinearity had to be used – temporal solitons can also be avoided in the focusing regime by considering the normal dispersion regime.

**Spatial vs spatio-temporal dynamics.** To clearly evidence the key role of the temporal degrees of freedom, we compare the ST simulation in Fig. 2, with the equivalent simulation in the purely spatial problem ($b_m(\omega, z) \rightarrow b_m^S(z)$, $\tilde{\beta}_m(\omega) \rightarrow \beta_m$, $\Gamma_m(\omega) \rightarrow 1$), where the spatial modal amplitudes are ruled by

$$i \partial_z b_m^S(z) = \beta_m b_m^S - \gamma \sum_{p,q,r} W_{mpqr} b_p^S b_q^{S*} b_r^S. \quad (7)$$

The evolution of the spatial mode distribution $N_m^S(z) = |b_m^S(z)|^2$ is then compared to the ST evolution $N_m^{ST}(z)$, considering the same initial condition and average power, $N_m^S(z=0) = N_m^{ST}(z=0)$. The comparison of the ST dynamics in Fig. 2, and the purely spatial dynamics in Fig. 3 is striking: The ST field exhibits a fast relaxation to equilibrium, whereas in the spatial case the thermalization process is frozen. This is confirmed by the evolution of the distance to equilibrium $D(z)$, which does not decrease (Fig. 3(b)), in contrast to the ST case in Fig. 2(b).

The evolution of the entropy also unveils the profound distinction between the ST and the pure spatial dynamics. In the ST case the irreversible process of thermalization is featured by a monotonic increase of entropy, as dictated by the H-theorem of entropy growth inherent to the kinetic Eq.(4), see Fig. 2(c). This is in contrast with the purely spatial simulation of Eq.(29), see Fig. 3(c): Here, the spatial dynamics evolves in a discrete wave turbulence regime governed by a formally reversible system of kinetic equations, which does not exhibit a H-theorem of entropy growth and that explains the frozen thermalization observed in Fig. 3 (see Fig. 5 in Appendix).

**Adiabatic cooling.** A well-known issue of classical field theories is the occurrence of UV divergences in the thermal equilibrium state, the so-called black-body catastrophe [52, 53, 68]. In our configuration, this issue is naturally tamed in the transverse direction by the finite number of modes of the waveguide, but gives rise to a rich physics in the temporal direction. Starting from the very non-thermal initial state with a short-tailed Gaussian distribution considered in our simulations, after a transient (typically $z \gtrsim 160L_{nl}$ in Fig. 2) the optical field approaches at each $z$ a local quasi-equilibrium state that closely approximates a RJ equilibrium distribution within a limited RJ window $[\tilde{\omega}_{c,-}^{loc}(z), \tilde{\omega}_{c,+}^{loc}(z)]$,
as illustrated in the mode-integrated spectrum shown in Fig.4(a). The accuracy of the local thermalization process is further evidenced by the remarkable agreement shown in Fig. 4(b) between the local RJ mode occupancies $N_m^{loc}(z)$, and the NEE simulation $N_m(z)$. For increasing $z$, the extension of the RJ window grows as shown in Fig. 4(c) and, correspondingly, the local thermodynamic parameters ($T^{loc}(z)$, $\lambda^{loc}(z)$, $\mu^{loc}(z)$) also display a marked $z$ dependence, as expected by imposing conservation of ($N_{\tilde{n}}$, $P_{\tilde{n}}$, $E_{\tilde{n}}$) in the presence of the $z$-dependent cut-offs $\tilde{\omega}_{c,\pm}^{loc}(z)$.

For instance, the temperature plotted in Fig. 4(d) displays a monotonic decrease due to the continuous transfer (at a constant energy $E_{\tilde{n}} = \text{const}$) of the incoherent beam fluctuations into the high-energy tails of the spectrum distribution. As a direct consequence of this conservative adiabatic cooling, a marked beam cleaning is visible in Fig. 4(b) as a transverse condensation effect with the population being macroscopically concentrated in the fundamental waveguide mode, $N_0 \gg N_{m \neq 0}$. Most interestingly, while in the simulations shown in this Letter the efficiency of the cooling process is limited by the numerical cut-off that one needs to impose to make the calculation feasible, in a real physical system the UV-divergent value of the energy stored in the tails of the classical distribution allows for an unbounded decrease of $T(z)$. Of course, this is only valid as long as a classical wave description is valid, that is before quantum statistical effects make the distribution to recover a quantum Bose-Einstein thermal law. Given the small value of typical nonlinear media, this is expected to occur at extremely long $z$ [49].

**Conclusions and outlook.** We have considered a most general description of broadband light propagation in nonlinear Kerr media to derive a wave turbulence formulation of spatio-temporal thermalization of a light beam propagating in a multi-mode waveguide. In contrast to the frozen spatial-only thermalization of monochromatic light due to the lack of efficient quasi-resonances among discrete transverse modes, the continuous nature of the temporal degrees leads to an efficient spatio-temporal thermalization. As a consequence of the blackbody catastrophe of classical fields, our route to ST thermalization unveils an intrinsic adiabatic cooling mechanism, whereby the field fluctuations are transferred to high-frequency components along the time dimension, so that the low-energy modes display a virtually unlimited spatial beam cleaning condensation. This adiabatic cooling is inherently conservative and therefore in contrast to conventional evaporative cooling techniques in Bose-Einstein condensates. Work in progress shows that the quasi-equilibrium process of adiabatic cooling can be described using the wave turbulence kinetic framework, providing valuable insights into the dynamics of nonequilibrium closed systems and their pathways to thermalization, see e.g. [75, 81]. From a broader perspective, our conservative wave-guided light configuration can eventually lead to a full 3D condensation in the quantum

regime, which would open novel avenues for ST beam cleaning and coherent light generation.

**Acknowledgments.** Funding was provided by Agence Nationale de la Recherche (Grants No. ANR-23-CE30-0021, No. ANR-21-ESRE-0040). Calculations were performed using HPC resources from Université Côte d'Azur's Center Azzura, and DNUM-CCUB (Université de Bourgogne Europe). I.C. acknowledges financial support by the PE0000023-NQSTI project by the Italian Ministry of University and Research, co-funded by the European Union - NextGeneration EU, and from Provincia Autonoma di Trento (PAT), partly via the Q@TN initiative. I.C.'s research was supported in part by the International Centre for Theoretical Sciences (ICTS) for participating in the program: Turbulence and Vortex dynamics in 2D quantum fluids (code: ICTS/QUFLU2024/2). A.P. thanks fruitful discussions with P. Béjot and B. Kibler.

# APPENDIX

## I. DERIVATION OF THE UPE (1-3)

We start from the Maxwell equations and consider the scalar approximation for the real-valued electric field $E(t, z, \mathbf{r})$. Considering the cubic (Kerr) nonlinearity, the Fourier transform of the field satisfies the nonlinear Helmholtz equation:

$$ \Delta \hat{E} + \frac{\omega^2}{c^2} \tilde{n}^2(\omega, \mathbf{r}) \hat{E} = -\frac{\omega^2}{c^2} \chi^{(3)} \hat{E}^3, \quad (8) $$

where $\Delta = \Delta_{\perp} + \partial_z^2$ is the three-dimensional Laplacian, and $c$ is the vacuum light speed. Here, $\tilde{n}(\omega, \mathbf{r})$ is the index of refraction of the waveguide at the frequency $\omega$. We denote $n(\omega) = \tilde{n}(\omega, \mathbf{0})$ the index of refraction at $\mathbf{r} = \mathbf{0}$, and $n_o = n(\omega_o)$, with $\omega_o$ the central frequency of the field. We may have $\tilde{n}^2(\omega, \mathbf{r}) = n^2(\omega)(1 - |\mathbf{r}|^2/a^2)$ for a parabolic (graded-index) fiber, or $\tilde{n}^2(\omega, \mathbf{r}) = n^2(\omega)$ for a homogeneous (step-index) fiber. We introduce the real-valued orthonormal basis $u_m(\mathbf{r})$, with eigenvalues $\beta_m$, which are solutions of

$$ -\frac{1}{2k_o} \Delta_{\perp} u_m + \frac{\omega_o}{2n_o c} (n_o^2 - n(\omega_o, \mathbf{r})^2) u_m = \beta_m u_m, $$

where $k_o = n_o \omega_o/c$. The field can be expanded on the basis $\hat{E}(\omega, z, \mathbf{r}) = \sum_m a_m(\omega, z) u_m(\mathbf{r})$, where $a_m(\omega, z) = \int u_m^*(\mathbf{r}) \hat{E}(\omega, z, \mathbf{r}) d\mathbf{r}$. Note that $a_m(-\omega, z) = a_m^*(\omega, z)$ because $E(t, z, \mathbf{r})$ is real-valued. By multiplying the Helmholtz equation (8) by $u_m$ and integrating in $\mathbf{r}$, we get

$$ \partial_z^2 a_m + (n_o^2 \omega^2 / c^2 - 2k_o \beta_m) a_m + \sum_p A_{mp}(\omega) a_p = - \frac{\chi^{(3)} \omega^2}{(2\pi)^2 c^2} \sum_{p,q,s} W_{mpqs} \int a_p(\omega_1) a_q(\omega_2) a_s(\omega_3) \delta_{\omega} d\omega_{123}, $$
with $W_{mpqs} = \int u_m(\mathbf{r})u_p(\mathbf{r})u_q(\mathbf{r})u_s(\mathbf{r})d\mathbf{r}$, $\hat{\delta}_\omega = \delta(\omega - \omega_1 - \omega_2 - \omega_3)$, $d\omega_{123} = d\omega_1 d\omega_2 d\omega_3$, and

$$A_{mp}(\omega) = \int Q(\omega, \mathbf{r}) u_m(\mathbf{r}) u_p(\mathbf{r}) d\mathbf{r},$$

$$Q(\omega, \mathbf{r}) = \frac{\omega^2}{c^2} (n^2(\omega, \mathbf{r}) - n_o^2) - \frac{\omega_o^2}{c^2} (n^2(\omega_o, \mathbf{r}) - n_o^2).$$

In general, the coefficients $a_m$ have rapid dominant phases of the form $\exp(ik_o z - i\beta_m z)$, so the terms $p \neq m$ in the sum over $p$ average out. The resulting equation can be written

$$\begin{aligned} \partial_z^2 a_m + (n_o^2 \omega^2 / c^2 - 2k_o \beta_m + \alpha_m(\omega)) a_m \\ = - \frac{\chi^{(3)} \omega^2}{(2\pi)^2 c^2} \sum_{p,q,s} W_{mpqs} \int a_p(\omega_1) a_q^*(\omega_2) a_s(\omega_3) \delta_\omega d\omega_{123}, \end{aligned} \quad (9)$$

with $\alpha_m(\omega) = A_{mm}(\omega)$. This equation is obtained without approximations for the case of interest of a homogeneous trapping potential (step-index fiber), since $A_{mp}(\omega) = \alpha(\omega)\delta_{mp}$ is diagonal, with $\alpha(\omega) = (\omega^2/c^2)(n(\omega)^2 - n_o^2)$. Note that we have changed $\omega_2 \to -\omega_2$, so that the Dirac function reads $\delta_\omega = \delta(\omega - \omega_1 + \omega_2 - \omega_3)$.

At variance with the multimode UPE in Ref. [67], here $\beta_m$ and $W_{mpqr}$ do not depend on the frequency $\omega$, because we have expanded the field on the eigenmodes at frequency $\omega_o$. Then we just need to compute $\beta_m$, $\alpha_m(\omega)$ and $W_{mpqs}$, and not a function indexed by $m, p, q, s$ and four frequencies. The latter would have been necessary if we had chosen to expand the field on $\omega$-dependent modes.

**The forward scattering approximation:** We first remark that the general solution of (9) without the nonlinear term has the form

$$a_m(\omega, z) = a_m^+(\omega) e^{i(k_o - \tilde{\beta}_m(\omega))z} + a_m^-(\omega) e^{-i(k_o - \tilde{\beta}_m(\omega))z},$$

with

$$\tilde{\beta}_m(\omega) = k_o - \sqrt{n_o^2 \omega^2 / c^2 - 2k_o \beta_m + \alpha_m(\omega)}. \quad (10)$$

Note that $\tilde{\beta}_m(\omega_o) \neq \beta_m$ in general. In the presence of the nonlinearity we introduce

$$a_m^+(\omega, z) = \frac{1}{2} \left( a_m(\omega, z) + \frac{\partial_z a_m(\omega, z)}{i(k_o - \tilde{\beta}_m(\omega))} \right) e^{-i(k_o - \tilde{\beta}_m(\omega))z},$$

$$a_m^-(\omega, z) = \frac{1}{2} \left( a_m(\omega, z) - \frac{\partial_z a_m(\omega, z)}{i(k_o - \tilde{\beta}_m(\omega))} \right) e^{i(k_o - \tilde{\beta}_m(\omega))z},$$

so that the solution of (9) has the form

$$a_m(\omega, z) = a_m^+(\omega, z) e^{i(k_o - \tilde{\beta}_m(\omega))z} + a_m^-(\omega, z) e^{-i(k_o - \tilde{\beta}_m(\omega))z},$$

and $a_m^+(\omega, z)$ and $a_m^-(\omega, z)$ satisfy the coupled first-order system

$$\partial_z a_m^+ = - \frac{i}{2(k_o - \tilde{\beta}_m(\omega))} F[\mathbf{a}]_m e^{-i(k_o - \tilde{\beta}_m(\omega))z}, \quad (11)$$

$$\partial_z a_m^- = + \frac{i}{2(k_o - \tilde{\beta}_m(\omega))} F[\mathbf{a}]_m e^{+i(k_o - \tilde{\beta}_m(\omega))z}. \quad (12)$$

Neglecting the backscattered components $a_m^-$, introducing the mode amplitudes defined by $b_m(\omega, z) = a_m(\omega, z)e^{-ik_o z}$, we obtain in a reference frame that propagates at the group velocity $v_g^{-1} = \partial_\omega k(\omega_o)$ the UPE (1):

$$i \partial_z b_m(\omega) = \tilde{B}_m(\omega) b_m - \gamma \Gamma_m(\omega) P_m(\mathbf{b}), \quad (13)$$

$$P_m(\mathbf{b}) = \sum_{pq s} W_{mpqs} \int b_p(\omega_1) b_q^*(\omega_2) b_s(\omega_3) \delta_\omega d\omega_{1-3}, \quad (14)$$

with $\tilde{B}_m(\omega) = \tilde{\beta}_m(\omega) + \frac{\omega}{v_g}$, $\Gamma_m(\omega) = \frac{\omega^2}{\omega_o^2} \frac{k_o}{k_o - \beta_m(\omega)}$ and $\gamma = \chi^{(3)}\omega_o / (8\pi^2 c)$.

The UPE (13-14) conserves the Hamiltonian $H_b = E_b + U_b$:

$$E_b = \sum_m \int \frac{\tilde{B}_m(\omega)}{\Gamma_m(\omega)} |b_m(\omega)|^2 d\omega, \quad (15)$$

$$U_b = -\frac{\gamma}{2} \sum_{mpqs} \int W_{mpqs} b_m^*(\omega_1) b_p(\omega_2) b_q^*(\omega_3) b_s(\omega_4) \delta_\omega d\omega_{1234} \quad (16)$$

with $\delta_\omega = \delta(\omega_1 + \omega_3 - \omega_2 - \omega_4)$. It also conserves the number of particles and the ‘temporal’ momentum

$$N_b = \sum_m \int \frac{|b_m(\omega)|^2}{\Gamma_m(\omega)} d\omega, \quad P_b = \sum_m \int \omega \frac{|b_m(\omega)|^2}{\Gamma_m(\omega)} d\omega. \quad (17)$$

The nonlinear length is defined by $L_{nl} = 1/(|\gamma_0|\bar{N}_b)$, where $\bar{N}_b = N_b/L$ is the optical power (averaged over the numerical time window, $L$), and $\gamma_0 = \gamma/A_{\text{eff}}^0$ the effective nonlinear coefficient (in $W^{-1}m^{-1}$), with $A_{\text{eff}}^0 = 1/W_{0000}$ the effective area of the fundamental mode.

**Application to a step-index waveguide:** The UPE (13,14) is general. It can be applied to a homogeneous (step-index) MMF, where $A_{mp}(\omega) = \alpha(\omega)\delta_{mp}$ with $\alpha(\omega) = (\omega^2/c^2)(n^2(\omega) - n^2(\omega_o))$. In this case $\tilde{\beta}_m(\omega)$ and $\Gamma_m(\omega)$ recover the expressions (2)-(3).

**The NLSE approximation:** Making use of the slowly-varying envelope approximation in the time domain and the paraxial approximation in the space domain, and using $k(\omega) = k(\omega_o) + \partial_\omega k(\omega_o)(\omega - \omega_o) + \frac{1}{2}\partial_\omega^2 k(\omega_o)(\omega - \omega_o)^2$, gives $\tilde{B}_m(\omega) = \beta_m - \kappa_2(\omega - \omega_o)^2$, because $\frac{1}{v_g} = \partial_\omega k(\omega_o)$ and $\kappa_2 = \frac{1}{2}\partial_\omega^2 k(\omega_o)$. The nonlinear coefficient reduces to a constant, $\Gamma_m(\omega) = 1$.

**The NEE approximation:** The main difference with respect to the NLSE limit is that the NEE preserves the frequency dependence of $k(\omega)$, but neglects the dispersion of the refractive index, i.e., we set $n(\omega) = n(\omega_o)$ in $k(\omega) = n(\omega_o)\omega/c$ [82], so that $\frac{2k_o\beta_m}{2k(\omega)} \approx \frac{\beta_m}{1 + \frac{\omega - \omega_o}{\omega_o}}$, and

$$\tilde{B}_m(\omega) = \frac{\beta_m}{1 + \frac{\omega - \omega_o}{\omega_o}} - \sum_{j \ge 2} \kappa_j (\omega - \omega_o)^j, \quad (18)$$
with $\kappa_j = \frac{1}{j} \partial_{\omega}^j k(\omega)|_{\omega_o}$. The first term in (18) couples spatial and temporal effects. The approximation $\beta_m(\omega) \simeq \beta_m(1 - \frac{\omega-\omega_o}{\omega_o}) - \sum_{j \ge 2} \kappa_j(\omega - \omega_o)^j$ is called partially-corrected NLS equation [82].

For the nonlinear coefficient $\Gamma_m(\omega)$, the NEE approximates $\sqrt{k(\omega)^2 - 2k_o\beta_m} \simeq k(\omega)$. Neglecting the dispersion of the refractive index, we get

$$ \Gamma_m(\omega) \approx \frac{\omega^2 k_o}{\omega_o^2 n(\omega_o)\omega/c} = 1 + \frac{\omega - \omega_o}{\omega_o}. $$

Note that the NEE can be written in the space-time do-

main. Considering the complex envelope field $\psi(\mathbf{r}, t, z)$, the NEE reads

$$ i \partial_z \psi(\mathbf{r}, t, z) = \frac{-1}{2k_o(1 + i\omega_o^{-1}\partial_t)} \nabla^2 \psi + \frac{1}{(1 + i\omega_o^{-1}\partial_t)} V(\mathbf{r})\psi + \sum_{j \ge 2} \kappa_j \partial_t^j \psi - \gamma(1 + i\omega_o^{-1}\partial_t) |\psi|^2 \psi. \quad (19) $$

By removing the operators $(1 + i\omega_o^{-1}\partial_t)$ and by retaining only second-order dispersion effects ($j = 2$), Eq. (19) recovers the NLSE approximation.

## II. DERIVATION OF THE WAVE TURBULENCE KINETIC EQ.(4-5)

We introduce a symmetric form of Eq. (13-14). The mode amplitudes $\tilde{b}_m(\omega, z) = b_m(\omega, z)/\sqrt{\Gamma_m(\omega)}$ satisfy

$$ i \partial_z \tilde{b}_m = \tilde{B}_m(\omega) \tilde{b}_m - \gamma \sum_{p,q,s} \iiint L_{mpqs}(\omega, \omega_1, \omega_2, \omega_3) \tilde{b}_p(\omega_1) \tilde{b}_q^*(\omega_2) \tilde{b}_s(\omega_3) \delta(\omega - \omega_1 + \omega_2 - \omega_3) d\omega_1 d\omega_2 d\omega_3, \quad (20) $$

with $L_{mpqs}(\omega, \omega_1, \omega_2, \omega_3) = W_{mpqs} \sqrt{\Gamma_m(\omega)\Gamma_p(\omega_1)\Gamma_q(\omega_2)\Gamma_s(\omega_3)}$, and recall that $\tilde{B}_m(\omega) = \tilde{\beta}_m(\omega) + \frac{\omega}{v_g}$. Equation (20) conserves the Hamiltonian $H_{\tilde{b}} = E_{\tilde{b}} + U_{\tilde{b}}$ with: $E_{\tilde{b}} = \sum_m \int \tilde{B}_m(\omega) |\tilde{b}_m(\omega)|^2 d\omega$, and

$$ U_{\tilde{b}} = -\frac{\gamma}{2} \sum_{mpqs} \iiint L_{mpqs}(\omega, \omega_1, \omega_2, \omega_3) \tilde{b}_m^*(\omega) \tilde{b}_p(\omega_1) \tilde{b}_q^*(\omega_2) \tilde{b}_s(\omega_3) \delta(\omega - \omega_1 + \omega_2 - \omega_3) d\omega d\omega_1 d\omega_2 d\omega_3 \quad (21) $$

and the number of particles $N_{\tilde{b}} = \sum_m \int |\tilde{b}_m(\omega)|^2 d\omega$, and ‘temporal’ momentum $P_{\tilde{b}} = \sum_m \int \omega |\tilde{b}_m(\omega)|^2 d\omega$.

We consider incoherent waves and we denote by $\langle \cdot \rangle$ the average over the realizations of the initial condition at $z = 0$. We have

$$ \begin{aligned} \partial_z \langle \tilde{b}_m(\omega_1) \tilde{b}_m^*(\omega_1') \rangle = i(\tilde{B}_m(\omega_1') - \tilde{B}_m(\omega_1)) \langle \tilde{b}_m(\omega_1) \tilde{b}_m^*(\omega_1') \rangle \\ + i\gamma \iiint \sum_{p,q,s} L_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4) \langle \tilde{b}_m^*(\omega_1') \tilde{b}_p(\omega_2) \tilde{b}_q^*(\omega_3) \tilde{b}_s(\omega_4) \rangle \delta(\omega_2 - \omega_3 + \omega_4 - \omega_1) d\omega_{2,3,4} \\ - i\gamma \iiint \sum_{p,q,s} L_{mpqs}^*(\omega_1, \omega_2, \omega_3, \omega_4) \langle \tilde{b}_m(\omega_1) \tilde{b}_p^*(\omega_2) \tilde{b}_q(\omega_3) \tilde{b}_s^*(\omega_4) \rangle \delta(\omega_2 - \omega_3 + \omega_4 - \omega_1') d\omega_{2,3,4}. \end{aligned} $$

and

$$ \begin{aligned} \partial_z \langle \tilde{b}_m^*(\omega_1) \tilde{b}_p(\omega_2) \tilde{b}_q^*(\omega_3) \tilde{b}_s(\omega_4) \rangle = i \Omega_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4) \langle \tilde{b}_m^*(\omega_1) \tilde{b}_p(\omega_2) \tilde{b}_q^*(\omega_3) \tilde{b}_s(\omega_4) \rangle \\ - i \gamma \sum_{p',q',s'} \iiint L_{mp'q's'}^*(\omega_1, \omega_1', \omega_2', \omega_3') \langle \tilde{b}_{p'}^*(\omega_1') \tilde{b}_{q'}(\omega_2') \tilde{b}_{s'}^*(\omega_3') \tilde{b}_p(\omega_2) \tilde{b}_q^*(\omega_3) \tilde{b}_s(\omega_4) \rangle \delta(\omega_1' - \omega_2' + \omega_3' - \omega_1) d\omega_{1,2,3} \\ + i \gamma \sum_{p',q',s'} \iiint L_{pp'q's'}(\omega_2, \omega_1', \omega_2', \omega_3') \langle \tilde{b}_m^*(\omega_1) \tilde{b}_{p'}(\omega_1') \tilde{b}_{q'}(\omega_2') \tilde{b}_{s'}^*(\omega_3') \tilde{b}_s(\omega_4) \rangle \delta(\omega_1' - \omega_2' + \omega_3' - \omega_2) d\omega_{1,2,3} \\ - i \gamma \sum_{p',q',s'} \iiint L_{qp'q's'}^*(\omega_3, \omega_1', \omega_2', \omega_3') \langle \tilde{b}_m^*(\omega_1) \tilde{b}_p(\omega_2) \tilde{b}_{p'}^*(\omega_1') \tilde{b}_{q'}(\omega_2') \tilde{b}_{s'}^*(\omega_3') \tilde{b}_s(\omega_4) \rangle \delta(\omega_1' - \omega_2' + \omega_3' - \omega_3) d\omega_{1,2,3} \\ + i \gamma \sum_{p',q',s'} \iiint L_{sp'q's'}(\omega_4, \omega_1', \omega_2', \omega_3') \langle \tilde{b}_m^*(\omega_1) \tilde{b}_p(\omega_2) \tilde{b}_q^*(\omega_3) \tilde{b}_{p'}(\omega_1') \tilde{b}_{q'}^*(\omega_2') \tilde{b}_{s'}(\omega_3') \rangle \delta(\omega_1' - \omega_2' + \omega_3' - \omega_4) d\omega_{1,2,3}, \end{aligned} $$

with $\Omega_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4) = \tilde{B}_m(\omega_1) - \tilde{B}_p(\omega_2) + \tilde{B}_q(\omega_3) - \tilde{B}_s(\omega_4)$. The equation for the evolution of the fourth-order moment (in the rhs of the above equation) will depend on the sixth-order moment. In this way, one obtains an infinite hierarchy of moment equations, in which the n-th order moment depends on the n + 2-order moment. The
hierarchy is closed in the weakly nonlinear regime because the field approaches Gaussian statistics. By virtue of the factorizability property of statistical Gaussian fields, we obtain

$$
\begin{align*}
\partial_z \langle \check{b}_m^*(\omega_1) \check{b}_p(\omega_2) \check{b}_q^*(\omega_3) \check{b}_s(\omega_4) \rangle &= i \Omega_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4) \langle \check{b}_m^*(\omega_1) \check{b}_p(\omega_2) \check{b}_q^*(\omega_3) \check{b}_s(\omega_4) \rangle \\
&\quad - 2i\gamma L_{mpqs}^*(\omega_1, \omega_2, \omega_3, \omega_4) \check{n}_m(\omega_1) \check{n}_p(\omega_2) \check{n}_q(\omega_3) \check{n}_s(\omega_4) \\
&\quad \times [\check{n}_m(\omega_1)^{-1} - \check{n}_p(\omega_2)^{-1} + \check{n}_q(\omega_3)^{-1} - \check{n}_s(\omega_4)^{-1}] \delta(\omega_3 - \omega_4 + \omega_1 - \omega_2) \\
&\quad - 2i\gamma \left[ \sum_{p'} \int L_{mpp'p'}^*(\omega_1, \omega_1, \omega, \omega) \check{n}_{p'}(\omega) d\omega \right] (\check{n}_p(\omega_1) - \check{n}_m(\omega_1)) \check{n}_q(\omega_3) \delta_{qs}^K \delta(\omega_2 - \omega_1) \delta(\omega_3 - \omega_4 + \omega_1 - \omega_2) \\
&\quad - 2i\gamma \left[ \sum_{p'} \int L_{mp'p's}^*(\omega_1, \omega, \omega, \omega_1) \check{n}_{p'}(\omega) d\omega \right] (\check{n}_s(\omega_1) - \check{n}_m(\omega_1)) \check{n}_p(\omega_2) \delta_{pq}^K \delta(\omega_4 - \omega_1) \delta(\omega_3 - \omega_4 + \omega_1 - \omega_2) \\
&\quad - 2i\gamma \left[ \sum_{p'} \int L_{p'p'qs}^*(\omega, \omega, \omega_3, \omega_3) \check{n}_{p'}(\omega) d\omega \right] (\check{n}_s(\omega_3) - \check{n}_q(\omega_3)) \check{n}_m(\omega_1) \delta_{pm}^K \delta(\omega_4 - \omega_3) \delta(\omega_3 - \omega_4 + \omega_1 - \omega_2) \\
&\quad - 2i\gamma \left[ \sum_{p'} \int L_{p'p'qp}^*(\omega, \omega, \omega_2, \omega_2) \check{n}_{p'}(\omega) d\omega \right] (\check{n}_p(\omega_2) - \check{n}_q(\omega_2)) \check{n}_m(\omega_1) \delta_{sm}^K \delta(\omega_3 - \omega_2) \delta(\omega_3 - \omega_4 + \omega_1 - \omega_2).
\end{align*}
$$

where $\langle \check{b}_m(\omega_1) \check{b}_{m'}^*(\omega_1') \rangle = \check{n}_m(\omega_1) \delta_{mm'}^K \delta(\omega_1 - \omega_1')$. By the presence of the factor $\delta(\omega_3 - \omega_4 + \omega_1 - \omega_2)$ in the rhs, we have

$$
\langle \check{b}_m^*(\omega_1) \check{b}_p(\omega_2) \check{b}_q^*(\omega_3) \check{b}_s(\omega_4) \rangle = \delta(\omega_3 - \omega_4 + \omega_1 - \omega_2) J_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4)
$$

and the spectrum $\check{n}_m(\omega_1)$ satisfies

$$
\partial_z \check{n}_m(\omega_1) = -2\gamma \sum_{p,q,s} \operatorname{Im} \left( \iiint L_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4) J_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4) \delta(\omega_1 + \omega_3 - \omega_2 - \omega_4) d\omega_{2,3,4} \right), \quad (22)
$$

$$
\begin{equation}
\begin{split}
\partial_z J_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4) = {}& i \Omega_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4) J_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_4) \\
& - 2i\gamma L_{mpqs}^*(\omega_1, \omega_2, \omega_3, \omega_4) M_{mpqs}[\mathbf{n}(z)](\omega_1, \omega_2, \omega_3, \omega_4) - 2i\gamma R_{mpqs}[\check{\mathbf{n}}(z)](\omega_1, \omega_2, \omega_3, \omega_4),
\end{split}
\tag{23}
\end{equation}
$$

where

$$
M_{mpqs}[\check{\mathbf{n}}](\omega_1, \omega_2, \omega_3, \omega_4) = \check{n}_m(\omega_1) \check{n}_p(\omega_2) \check{n}_q(\omega_3) \check{n}_s(\omega_4) (\check{n}_m^{-1}(\omega_1) + \check{n}_q^{-1}(\omega_3) - \check{n}_p^{-1}(\omega_2) - \check{n}_s(\omega_4)^{-1}), \quad (24)
$$

$$
\begin{align*}
R_{mpqs}[\check{\mathbf{n}}](\omega_1, \omega_2, \omega_3, \omega_4) &= U_{mp}^*[\check{\mathbf{n}}](\omega_1) (\check{n}_p(\omega_2) - \check{n}_m(\omega_1)) \check{n}_q(\omega_3) \delta_{qs}^K \delta(\omega_4 - \omega_3) \\
&\quad + U_{ms}^*[\check{\mathbf{n}}](\omega_1) (\check{n}_s(\omega_4) - \check{n}_m(\omega_1)) \check{n}_p(\omega_2) \delta_{pq}^K \delta(\omega_3 - \omega_2) \\
&\quad + U_{qs}^*[\check{\mathbf{n}}](\omega_3) (\check{n}_s(\omega_4) - \check{n}_q(\omega_3)) \check{n}_m(\omega_1) \delta_{pm}^K \delta(\omega_2 - \omega_1) \\
&\quad + U_{qp}^*[\check{\mathbf{n}}](\omega_3) (\check{n}_p(\omega_2) - \check{n}_q(\omega_3)) \check{n}_m(\omega_1) \delta_{sm}^K \delta(\omega_4 - \omega_1),
\end{align*}
$$

$$
\mathcal{U}_{qp}[\check{\mathbf{n}}](\omega) = \sum_{p'} \int L_{qpp'p'}(\omega, \omega, \omega', \omega') \check{n}_{p'}(\omega') d\omega'. \quad (25)
$$

After simplification, we finally obtain

$$
\partial_z \check{n}_m(\omega_1) = -2\gamma \operatorname{Im} \left( \sum_{p,q,s} \iint J_{mpqs}^{(1)}(\omega_1, \omega_2, \omega_3) d\omega_{2,3} + 2 \sum_p J_{mp}^{(2)}(\omega_1) \right), \quad (26)
$$

with

$$
\begin{equation}
\begin{split}
\partial_z J_{mpqs}^{(1)}(\omega_1, \omega_2, \omega_3) = {}& i \Omega_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_1 + \omega_3 - \omega_2) J_{mpqs}^{(1)}(\omega_1, \omega_2, \omega_3) \\
& - 2i\gamma |L_{mpqs}(\omega_1, \omega_2, \omega_3, \omega_1 + \omega_3 - \omega_2)|^2 M_{mpqs}[\check{\mathbf{n}}(z)](\omega_1, \omega_2, \omega_3, \omega_1 + \omega_3 - \omega_2),
\end{split}
\tag{27}
\end{equation}
$$

and

$$
\partial_z J_{mp}^{(2)}(\omega_1) = i(\tilde{B}_m(\omega_1) - \tilde{B}_p(\omega_1)) J_{mp}^{(2)}(\omega_1) - 2i\gamma |\mathcal{U}_{mp}[\tilde{\mathbf{n}}(z)](\omega_1)|^2 (\tilde{n}_p(\omega_1) - \tilde{n}_m(\omega_1)). \quad (28)
$$

Eqs. (26-28) are the equations driving the evolution of the spectrum. Making use of the wave turbulence theory [53, 70], and neglecting degenerate modes, one obtains the kinetic Eq.(4) for the evolution of the modal spectra $\check{n}_m(\omega, z)$.
### III. PURE SPATIAL CASE: COUPLED SECOND- AND FOURTH-ORDER MOMENTS EQUATIONS

**Moments equations:** The starting point is the pure spatial model Eq.(7):

$$i \partial_z b_m^S(z) = \beta_m b_m^S - \gamma \sum_{p,q,r} W_{mpqr} b_p^S b_q^{S*} b_r^S. \qquad (29)$$

where linear and nonlinear dispersion relations no longer depend on the frequency: $\tilde{\beta}_m(\omega) \to \beta_m$, $\gamma\Gamma_m(\omega) \to \gamma$, and then $b_m(\omega, z) \to b_m^S(z)$. Here, we show that when nonresonant interactions dominate the dynamics, the wave turbulence kinetic equations do not exhibit a H-theorem of entropy growth.

Considering the pure spatial dynamics of Eq.(29), the coupled Eqs.(22-23) for the second and fourth order moments

$$n_m^S(z) = \langle b_m^S b_m^{S*} \rangle, \quad J_{mpqs}^S(z) = \langle b_m^{S*} b_p^S b_q^{S*} b_s^S \rangle,$$

reduce to

$$\partial_z n_m^S = +\gamma i \sum_{p,q,s} W_{mpqs} J_{mpqs}^S - \gamma i \sum_{p,q,s} W_{mpqs}^{*} J_{mpqs}^{S*} \qquad (30)$$

$$\partial_z J_{mpqs}^S = i \Delta \beta_{mpqs} J_{mpqs}^S - 2i\gamma W_{mpqs}^{*} M_{mpqs}(n(z)) - 2i\gamma R_{mpqs}(n(z)) \qquad (31)$$

where $\Delta\beta_{mpqs} = \beta_m + \beta_q - \beta_p - \beta_s$, and

$$M_{mpqs}(\mathbf{n}) = n_p^S n_q^S n_s^S + n_m^S n_p^S n_s^S - n_m^S n_q^S n_s^S - n_m^S n_p^S n_q^S,$$

$$R_{mpqs}(\mathbf{n}) = \delta_{q,s}^K U_{pm}(\mathbf{n})(n_p^S - n_m^S)n_q^S + \delta_{q,p}^K U_{sm}(\mathbf{n})(n_s^S - n_m^S)n_p^S \delta_{m,s}^K U_{pq}(\mathbf{n})(n_p^S - n_q^S)n_s^S + \delta_{m,p}^K U_{sq}(\mathbf{n})(n_s^S - n_q^S)n_m^S,$$

$$U_{pq}(\mathbf{n}) = \sum_{s} W_{pqss} n_s^S = \int u_p^*(\mathbf{r}) u_q(\mathbf{r}) \sum_{s} n_s^S |u_s(\mathbf{r})|^2 d\mathbf{r}.$$

Eq.(30-31) conserve the power (number of particles) and the total energy:

$$N = \sum_{m} n_{m}^{S}(z), \quad H = \sum_{m} \beta_{m} n_{m}^{S}(z) - \frac{\gamma}{4} \sum_{mpqs} W_{mpqs} J_{mpqs}^{S}(z) + W_{mpqs}^{*} J_{mpqs}^{S*}(z).$$

In addition, Eq.(30-31) are formally reversible, i.e., they are invariant under the change of variable

$$z \to -z, \quad J_{mpqs}^S \to J_{mpqs}^{S*}, \quad W_{mpqs} \to W_{mpqs}^{*}.$$

Consequently Eq.(30-31) do not exhibit a H-theorem of entropy growth.

**Numerical simulations of Eqs.(30-31) explaining the frozen thermalization of the pure spatial dynamics:** We have performed numerical simulations of the coupled second-order and fourth-order moments Eqs.(30-31), starting from the same initial condition as in Fig. 2. We compare the numerical results to those of the simulations of the purely spatial model Eq.(29) reported in Fig. 3, where an average over 21 realizations was taken. The results are reported in Fig. 5. We observe a quantitative agreement without using adjustable parameters. The kinetic Eqs.(30-31) then explain the frozen process of thermalization observed in Fig. 3.

**Discussion:** In order to derive the classical irreversible kinetic equation for the second-order moment, one needs to take the continuous limit of the discrete sums over the modes in Eq.(30-31). Before taking any limit, we note that the sums in Eq.(30-31) involve three different types of resonances:

(i) The exact resonances correspond to combinations of the uples $\{m, p, q, s\}$ such that $\Delta\beta_{mpqs} = 0$. These include in particular all trivial resonances involving only two spatial modes, i.e., $(m = p, q = s)$, or $(m = s, q = p)$. The formal solution of Eq.(31) for $J_{mpqs}$ can be substituted in Eq.(30), which gives for $\Delta\omega_{mpqs} = 0$:

$$\partial_z^2 n_m^S = 4\gamma^2 \sum_{p,q,s} \delta^K(\Delta\beta_{mpqs}) |W_{mpqs}|^2 M_{mpqs}(n) + 4\gamma^2 \sum_{p,q,s} \delta^K(\Delta\beta_{mpqs}) \Re[W_{mpqs} R_{mpqs}(n(z))]. \qquad (32)$$

This equation is formally reversible (second-order with respect to the 'time' z-variable). Exact resonances occur in the case where the trapping potential $V(r)$ exhibits a parabolic shape, due to the regular spacing of the eigenvalues. In this case, the simulations evidence a quasi-reversible exchange of power among the modes, which tends to freeze
thermalization process. This oscillatory behavior is related to the existence of Fermi-Pasta-Ulam recurrences, as discussed in Ref.[83] in the weakly nonlinear regime of the 2D NLS equation with a parabolic trapping potential. Note that in this latter case, RJ thermalization can be restored by introducing a weak disorder (random mode coupling) that breaks the coherent modal phase dynamics [14, 84].

(ii) The non-resonant terms correspond to combinations of the uples {$m, p, q, s$} such that $|\Delta\beta_{mpqs}|L_{nl} \gg 1$. These terms are characterized by a rapid rotating phase of $J_{mpqs}(z)$ that averages to zero the evolutions of $n_m(z)$.

(iii) The quasi-resonant terms are the uples {$m, p, q, s$} such that $|\Delta\beta_{mpqs}|L_{nl} \lesssim 1$.

If quasi-resonances dominate over exact resonances, then the dynamics is irreversible. This happens in the continuous limit, because there are many more quasi-resonances than exact resonances. In this limit, one recovers the classic kinetic equation that exhibits a H-theorem of entropy growth describing RJ thermalization. However, the continuous limit is not justified when one considers usual optical experiments with highly multimode step-index fibers.

<figure><img src="image_6.png" /><figcaption>FIG. 5: Spatial case: Frozen thermalization. The blue lines report the numerical simulation of the kinetic Eqs.(30-31) showing the spectrum $n_m^S(z)$ at $z = 600L_{nl}$ (a), and corresponding evolutions of the distance to the RJ equilibrium $D(z)$ (b), and the entropy $S(z)$ (c). The orange lines in (a)-(b)-(c) report the numerical simulations of Eq.(29) governing the spatial modal amplitudes $b_m^S(z)$: Because of the large fluctuations, an average over 21 realizations has been taken by starting from the same initial spectrum (solid black line), with different realizations of the random phases. A good agreement between the kinetic Eqs.(30-31) and the spatial model Eq.(29) is obtained without using adjustable parameters. The dashed black line in (a) reports the expected RJ equilibrium spectrum. The kinetic Eq.(30-31) then explains the frozen thermalization of the purely spatial dynamics $b_m^S(z)$ [Eq.(29)] discussed in Fig. 3, as confirmed by the evolutions of the distance $D(z)$ and the entropy $S(z)$ in panels (b)-(c).</figcaption></figure>

## IV. NUMERICAL METHODS

The ST simulations of the UPE (1) have been performed in the NEE approximation. Considering the time-frequency window used in the simulations ($\tilde{\omega}_c\tau_0 = 40$, see Eq.(36)), there is no appreciable difference between the UPE and the NEE approximation for the linear and nonlinear dispersion relations, $\tilde{\beta}_p(\omega)$ and $\Gamma_p(\omega)$, respectively. The NEE equation is solved using a pseudo-spectral split-step method, with a Fourier truncated spectrum defined through a Galerkin truncation [54]. In order to accurately conserve the momentum $P_b$, we have implemented

a dealiasing numerical procedure. In this way, $P_b$ and $N_b$ in Eqs.(17) are conserved at $10^{-7}$, and the Hamiltonian $H_b$ in Eqs.(15-16) at $10^{-4}$, throughout the simulation reported in Fig. 2. We considered a step index waveguide that guides 26 modes with an ellipticity that removes the mode degeneracies, in order to be consistent with the derivation of the wave turbulence kinetic equation. The propagation constants lie within the interval $\beta_0 = 25/L_{nl}$ and $\beta_{25} \simeq 500/L_{nl}$, so that the spatial mode dynamics evolves in the weakly nonlinear regime, $\beta_p L_{nl} \gg 1$. In Fig. 2 we considered the anomalous dispersion regime ($\kappa_2 < 0$) with a defocusing nonlinearity ($\gamma < 0$), to avoid the formation of temporal solitons, whose presence slows, or even freezes, the thermalization process. The simulations are realized in dimensionless units, and corresponding typical parameters can be $\kappa_2 = -0.1\text{ps}^2/\text{m}$, $\kappa_3 = 1.8 \times 10^{-4}\text{ps}^3/\text{m}$, $L_{nl} = 0.4\text{m}$. The presence of temporal solitons can also be avoided in the focusing regime ($\gamma > 0$) by considering the normal dispersion regime ($\kappa_2 > 0$).

### A. Global RJ equilibrium

In this section, we discuss the computation of the global RJ equilibrium distribution throughout the time-frequency window available in the simulation, as reported in Fig. 2. The case of local RJ equilibrium where the equilibrium distribution is computed over a reduced time-frequency window (Fig. 4) will be discussed in the next subsection.

The three parameters ($T, \lambda, \mu$) involved in the RJ distribution (6) are determined by the three conserved quantities ($E_b, P_b, N_b$):

$$E_b = TM + \lambda P_b + \mu N_b \quad (33)$$

$$\frac{P_b}{N_b} = \frac{\sum_{m,j} \frac{\tilde{\omega}_j}{\beta_m(1+\tilde{\omega}_j/\omega_o)^{-1} - \sum_{k \ge 2} \kappa_k \tilde{\omega}_j^k - \lambda \tilde{\omega}_j - \mu}}{\sum_{m,j} \frac{1}{\beta_m(1+\tilde{\omega}_j/\omega_o)^{-1} - \sum_{k \ge 2} \kappa_k \tilde{\omega}_j^k - \lambda \tilde{\omega}_j - \mu}} \quad (34)$$

$$\frac{E_b}{N_b} = \frac{\sum_{m,j} \frac{\beta_m(1+\tilde{\omega}_j/\omega_o)^{-1} - \sum_{k \ge 2} \kappa_k \tilde{\omega}_j^k}{\beta_m(1+\tilde{\omega}_j/\omega_o)^{-1} - \sum_{k \ge 2} \kappa_k \tilde{\omega}_j^k - \lambda \tilde{\omega}_j - \mu}}{\sum_{m,j} \frac{1}{\beta_m(1+\tilde{\omega}_j/\omega_o)^{-1} - \sum_{k \ge 2} \kappa_k \tilde{\omega}_j^k - \lambda \tilde{\omega}_j - \mu}} \quad (35)$$

where $M = qp$ is the total number of modes in the nu-
mericas (q spatial modes, p temporal modes). We recall that $\tilde{\omega} = \omega - \omega_o \ll \omega_o$ is the small frequency offset with respect to the central frequency $\omega_o$. The mode time-frequencies $\tilde{\omega}_j = j2\pi/L$ ($j = -p/2, -1, ..., p/2$, with p modes) result from the temporal grid used in the simulations (with periodic boundary conditions), where L is the time window. Accordingly, the frequency cut-off for the positive and negative frequencies are (approximately) the same, and they will be denoted hereafter by $\tilde{\omega}_c$:

$$ \tilde{\omega}_c = \tilde{\omega}_{c,+} = \pi p / L \approx -\tilde{\omega}_{c,-}. \quad (36) $$

The two parameters ($\lambda, \mu$) solutions of (34,35) are determined by an optimization algorithm. Then T is obtained from (33). Note that, by increasing the number of modes p (while keeping fixed the frequency cut-off $\tilde{\omega}_c = \tilde{\omega}_p$), the discrete sums over the temporal modes $\tilde{\omega}_j$ can be converted to continuous integrals, and the parameters ($T, \lambda, \mu$) solution of Eqs.(33-35) converge to well-defined values, that is, the RJ equilibrium distribution (6) does not depend on the numerical discretization of the temporal grid. The calculations reveal that the number of modes used in the simulations ($p = 256$) is sufficient to reach the continuous limit (with $q = 26$ spatial modes). The limited number of temporal modes results from the significant CPU time required for the long propagation lengths $L \sim 2 \times 10^6 \beta_{\text{max}}^{-1}$ necessary to reach thermal equilibrium ($\beta_{\text{max}}^{-1}$ being the smallest propagation length scale).

## B. Local RJ equilibrium

As discussed through Fig. 4, the optical field follows a local RJ equilibrium state during the propagation, where the parameters ($T^{\text{loc}}(z)$, $\lambda^{\text{loc}}(z)$, $\mu^{\text{loc}}(z)$) are computed over a reduced time frequency window $\tilde{\omega} \in [\tilde{\omega}_{c,-}^{\text{loc}}, \tilde{\omega}_{c,+}^{\text{loc}}]$, which is characterized by two local frequency cut-off $\tilde{\omega}_{c,\pm}^{\text{loc}}$, for the positive and negative frequencies:

$$ -\tilde{\omega}_c \le \tilde{\omega}_{c,-}^{\text{loc}} < 0 < \tilde{\omega}_{c,+}^{\text{loc}} \le \tilde{\omega}_c. \quad (37) $$

Note that $\tilde{\omega}_{c,+}^{\text{loc}}(z) \neq -\tilde{\omega}_{c,-}^{\text{loc}}(z)$ because of the intrinsic asymmetry of the spectral broadening process associated with the NEE model – the spectrum is symmetric for the more simple NLSE model. The fact that the field follows a local RJ equilibrium state during propagation is supported by the following numerical analysis: (i) At some propagation length $z_0$, a local average of the optical field is computed over $z \in [z_0 - \Delta z, z_0 + \Delta z]$, with $\Delta z = 3L_{nl}$ (10 realizations). (ii) At $z = z_0$, we compute the two local frequency cut-off, $\tilde{\omega}_{c,\pm}^{\text{loc}}(z_0)$, that minimize the distance $\mathcal{D}^{\text{loc}}(z_0) = \sum_p |N_p(z_0) - N_p^{\text{RJ,loc}}| / \sum_p (N_p(z_0) + N_p^{\text{RJ,loc}})$, between the averaged optical field and the local RJ equilibrium, see Fig. 6(a). This provides the local thermodynamic parameters: ($T^{\text{loc}}(z_0)$, $\lambda^{\text{loc}}(z_0)$, $\mu^{\text{loc}}(z_0)$) that characterize the local RJ equilibrium at $z = z_0$ over the reduced frequency window $[\tilde{\omega}_{c,-}^{\text{loc}}, \tilde{\omega}_{c,+}^{\text{loc}}]$, see Figs. 6(e)-(f) and Fig. 4(d).

<figure><img src="image_7.png" /><figcaption>FIG. 6: Local-equilibrium route to ST thermalization. (a) Local distance $\mathcal{D}^{\text{loc}}(z)$ computed over the reduced frequency interval $\tilde{\omega} \in [\tilde{\omega}_{c,-}^{\text{loc}}, \tilde{\omega}_{c,+}^{\text{loc}}]$ (blue line), and global distance $\mathcal{D}(z)$ computed over the whole frequency window $[-\tilde{\omega}_c, \tilde{\omega}_c]$ (orange line). Panels (b-d) report the mode-integrated temporal spectrum of the field from the NEE simulation at different propagation lengths (blue line), and the corresponding local RJ equilibrium distribution computed over the reduced frequency window $[\tilde{\omega}_{c,-}^{\text{loc}}, \tilde{\omega}_{c,+}^{\text{loc}}]$ (orange). Panels (e-f) report the corresponding evolution through propagation of the local thermodynamic parameters $\lambda^{\text{loc}}(z)$ and $\mu^{\text{loc}}(z)$ that characterize the local RJ equilibrium – the corresponding evolution of the temperature $T^{\text{loc}}(z)$ is reported in Fig. 4(d). (g) Evolution during propagation of the normalized second-order moment $\Omega_m(z)$ of the time-spectrum of the field for different modes $m$, which are retrieved from the NEE simulation (solid lines), and from the local RJ equilibrium spectrum (dashed lines). Note that the fundamental mode $m = 0$ shows spectral narrowing, signaling the onset of a ST beam cleaning process.</figcaption></figure>

We report in Figs. 6(b)-(d) the mode-integrated field temporal spectrum $\sum_m |b_m(\omega)|^2$ of the NEE simulation (Fig. 2), and the corresponding local RJ equilibrium spectrum, at different propagation lengths. At variance with the negative frequency cutoff $\tilde{\omega}_{c,-}^{\text{loc}}$, the positive frequency cut-off increases rapidly and reaches the cut-off frequency of the spectral grid considered in the simulation, that is $\tilde{\omega}_{c,+}^{\text{loc}} \approx \tilde{\omega}_c = 40/\tau_0$, see Fig. 4(a). The fact that the optical field follows a quasi-equilibrium path through propagation is reflected by the small value of the distance to the local RJ equilibrium, $\mathcal{D}^{\text{loc}}(z_0) < 0.05$ in Fig. 6(a). For clarity, we have also plotted in Fig. 6(a) the global distance $\mathcal{D}(z)$ to the global RJ reported in Fig. 2(b), which is computed with the frequency cut-off.
off related to the numerical spectral grid, $\tilde{\omega}_c = 40/\tau_0$. We stress in Fig. 4(b) the remarkable agreement between the spatial modal population computed from the local RJ equilibrium $N_m^{loc}(z)$ and the actual modal population $N_m(z)$ in the simulation. The fitting procedure is based on minimizing the distance $D^{loc}(z)$, which means that nothing in the procedure explicitly enforces the agreement with the evolution of individual modes, $N_m^{loc}(z)$ and $N_m^{ST}(z)$, shown in Fig. 4(b).

Finally, we report in Fig. 6(g) the evolution during propagation of the normalized second-order moment of the spectrum:

$$ \Omega_m(z) = \int (\omega - \omega_o)^2 |b_m(\omega, z)|^2 d\omega / \int |b_m(\omega, z)|^2 d\omega. $$

The solid lines report the evolution of $\Omega_m(z)$ recovered from the NEE simulation (Fig. 2), the dashed lines the evolution retrieved from the local RJ equilibrium spectra. It is interesting to note that the fundamental mode is the unique mode whose second-order moment decreases during propagation, indicating an incipient process of spatio-temporal beam cleaning.

[1] I. Carusotto, Quantum Fluids of Light, in Encyclopedia of Condensed Matter Physics, 2nd edition. Available at arXiv:2211.10980.

[2] J. Bloch, I. Carusotto, M. Wouters, Nonequilibrium Bose-Einstein condensation in photonic systems, Nature Phys. Rev. 4, 470 (2022).

[3] G. Situ, J.W. Fleischer, Dynamics of the Berezinskii-Kosterlitz-Thouless transition in a photon fluid, Nature Photonics 14, 517 (2020).

[4] Q. Fontaine, T. Bienaimé, S. Pigeon, E. Giacobino, A. Bramati, and Q. Glorieux, Observation of the Bogoliubov Dispersion in a Fluid of Light, Phys. Rev. Lett. 121, 183604 (2018).

[5] C. Michel, O. Boughdad, M. Albert, P.-E. Larré, and M. Bellec, Superfluid motion and drag-force cancellation in a fluid of light, Nature Comm. 9, 2108 (2018).

[6] M. Abuzarli, N. Cherroret, T. Bienaimé, Q. Glorieux, Nonequilibrium Prethermal States in a Two-Dimensional Photon Fluid, Phys. Rev. Lett. 129, 100602 (2022).

[7] R. Bekenstein, R. Schley, M. Mutzafi, C. Rotschild, and M. Segev, Optical simulations of gravitational effects in the Newton-Schrödinger system, Nature Phys. 11, 872 (2015).

[8] M.C. Braidotti, R. Prizia, C. Maitland, F. Marino, A. Prain, I. Starshynov, N. Westerberg, E.M. Wright, and D. Faccio, Measurement of Penrose superradiance in a photon superfluid, Phys. Rev. Lett. 128, 013901 (2022).

[9] G. Marcucci and C. Conti, Simulating general relativity and non-commutative geometry by non-paraxial quantum fluids, New Journal of Physics 21, 123038 (2019).

[10] L.G. Wright, Z. Liu, D.A. Nolan, M.-J. Li, D.N. Christodoulides, F.W. Wise, Self-organized instability in graded-index multimode fibres, Nature Photonics 10, 771 (2016).

[11] K. Krupa, A. Tonello, B.M. Shalaby, M. Fabert, A. Barthélémy, G. Millot, S. Wabnitz, V. Couderc, Spatial beam self-cleaning in multimode fibres, Nature Photonics 11, 237 (2017).

[12] M. Ferraro, F. Mangini, M. Zitelli, S. Wabnitz, On spatial beam self-cleaning from the perspective of optical wave thermalization in multimode graded-index fibers, Advances in Physics: X 8, 2228018 (2023).

[13] F.O. Wu, A.U. Hassan, D.N. Christodoulides, Thermo-

dynamic theory of highly multimoded nonlinear optical systems, Nature Photonics 13, 776 (2019).

[14] A. Fusaro, J. Garnier, K. Krupa, G. Millot, A. Picozzi, Dramatic acceleration of wave condensation mediated by disorder in multimode fibers, Phys. Rev. Lett. 122, 123902 (2019).

[15] J. Garnier, A. Fusaro, K. Baudin, C. Michel, K. Krupa, G. Millot, A. Picozzi, Wave condensation with weak disorder versus beam self-cleaning in multimode fibers, Phys. Rev. A 100, 053835 (2019).

[16] E.V. Podivilov, F. Mangini, O.S. Sidelnikov, M. Ferraro, M. Gervaziev, D.S. Kharenko, M. Zitelli, M.P. Fedoruk, S.A. Babin, S. Wabnitz, Thermalization of orbital angular momentum beams in multimode optical fibers, Phys. Rev. Lett. 128, 243901 (2022).

[17] A. Ramos, L. Fernández-Alcázar, T. Kottos, B. Shapiro, Optical Phase Transitions in Photonic Networks: a Spin-System Formulation, Phys. Rev. X 10, 031024 (2020).

[18] F.O. Wu, Q. Zhong, H. Ren, P.S. Jung, K.G. Makris, D.N. Christodoulides, Thermalization of Light's orbital angular momentum in nonlinear multimode waveguide systems, Phys. Rev. Lett. 128, 123901 (2022).

[19] Q. Zhong, F.O. Wu, A.U. Hassan, R. El-Ganainy, D.N. Christodoulides, Universality of light thermalization in multimoded nonlinear optical systems, Nature Comm. 14, 370 (2023).

[20] K. Baudin, A. Fusaro, K. Krupa, J. Garnier, S. Rica, G. Millot, A. Picozzi, Classical Rayleigh-Jeans condensation of light waves: Observation and thermodynamic characterization, Phys. Rev. Lett. 125, 244101 (2020).

[21] H. Pourbeyram, P. Sidorenko, F. Wu, L. Wright, D. Christodoulides, F. Wise, Direct observations of thermalization to a Rayleigh-Jeans distribution in multimode optical fibres, Nature Physics 18, 685 (2022).

[22] F. Mangini, M. Gervaziev, M. Ferraro, D. S. Kharenko, M. Zitelli, Y. Sun, V. Couderc, E.V. Podivilov, S.A. Babin, and S. Wabnitz, Statistical mechanics of beam self-cleaning in GRIN multimode optical fibers, Opt. Exp. 30, 10850 (2022).

[23] M. Zitelli, F. Mangini, S. Wabnitz, Statistics of modal condensation in nonlinear multimode fibers, Nature Comm. 15, 1149 (2024).

[24] A.L. Muniz, F. Wu, P. Jung, M. Khajavikhan, D.
Christodoulides, U. Peschel, Observation of photon-photon thermodynamic processes under negative optical temperature conditions, Science **379**, 1019 (2023).

[25] K. Baudin, J. Garnier, A. Fusaro, N. Berti, C. Michel, K. Krupa, G. Millot, A. Picozzi, Observation of light thermalization to negative temperature Rayleigh-Jeans equilibrium states in multimode optical fibers, Phys. Rev. Lett. **130**, 063801 (2023).

[26] M. Ferraro, F. Mangini, K. Stefánska, W.A. Gemechu, F. Frezza, V. Couderc, M. Gervaziev, D. Kharenko, S. Babin, and S. Wabnitz, Negative absolute temperature attractor in a dense photon gas, arXiv:2505.21163

[27] M. Ferraro, F. Mangini, F.O. Wu, M. Zitelli, D.N. Christodoulides, S. Wabnitz, Calorimetry of photon gases in nonlinear multimode optical fibers, Phys. Rev. X **14**, 021020 (2024).

[28] M.S. Kirsch, G.G. Pyrialakos, R. Altenkirch, M.A. Selim, J. Beck, T.A.W. Wolterink, H. Ren, P.S. Jung, M. Kha-javikhan, A. Szameit, M. Heinrich, D.N. Christodoulides, Observation of Joule-Thomson photon-gas expansion, Nat. Phys. **21**, 214-220 (2025).

[29] N. Berti, K. Baudin, A. Fusaro, G. Millot, A. Picozzi, J. Garnier, Interplay of thermalization and strong disorder: Wave turbulence theory, numerical simulations, and experiments in multimode optical fibers, Phys. Rev. Lett. **129**, 063901 (2022).

[30] M. Lian, Y. Geng, Y.-J. Chen, Y. Chen, J.-T. Lü, Coupled Thermal and Power Transport of Optical Waveguide Arrays: Photonic Wiedemann-Franz Law and Rectification Effect, Phys. Rev. Lett. **133**, 116303 (2024).

[31] A. Kurnosov, L. J. Fernández-Alcázar, A. Ramos, B. Shapiro, T. Kottos, Optical Kinetic Theory of Nonlinear Multimode Photonic Networks, Phys. Rev. Lett. **132**, 193802 (2024).

[32] Y. Liu, J. Lu, Z. Xiong, F.O. Wu, D. Christodoulides, Y. Chen, J.H. Jiang, Nonequilibrium transport and the fluctuation theorem in the thermodynamic behaviors of nonlinear photonic systems, Phys. Rev. Res. **7**, 013084 (2025).

[33] P.-E. Larré, I. Carusotto, Propagation of a quantum fluid of light in a cavityless nonlinear optical medium: General theory and response to quantum quenches, Physical Review A **92**, 043802 (2015).

[34] J. Steinhauer, M. Abuzarli, T. Aladjidi, T. Bienaimé, C. Piekarski, W. Liu, E. Giacobino, A. Bramati, Q. Giorrieux, Analogue cosmological particle creation in an ultracold quantum fluid of light, Nature Comm. **13**, 2890 (2022).

[35] A. Couairon, D. Faccio, and P. Di Trapani, Conical Waves, Filaments and Nonlinear Filamentation Optics (Aracne, Rome, 2007).

[36] B.A. Malomed, D. Mihalache, F. Wise, L. Torner, Spatiotemporal optical solitons, J. Opt. B: Quantum Sem class. Opt. **7**, R53 (2005).

[37] S. Chen, J.M. Dudley, Spatiotemporal Nonlinear Optical Self-Similarity in Three Dimensions, Phys. Rev. Lett. **102**, 233903 (2009).

[38] L.G. Wright, F.O. Wu, D.N. Christodoulides, F.W. Wise, Physics of highly multimode nonlinear optical systems, Nature Physics **18**, 1018 (2022).

[39] G. Agrawal, *Nonlinear Fiber Optics* (Academic, New York, 6th ed., 2019).

[40] Y. Sun, P. Parra-Rivas, G.P. Agrawal, T. Hansson, C. Antonelli, A. Mecozzi, F. Mangini, S. Wabnitz, Multi-

mode solitons in optical fibers: a review, Photon. Research **12**, 2581 (2024).

[41] Spatiotemporal mode-locking in fiber lasers represents another significant area of research. However, we do not address this topic as it pertains to forced-dissipative systems, in contrast to conservative (Hamiltonian) systems discussed in this work.

[42] B. Kibler, P. Béjot, Discretized Conical Waves in Multimode Optical Fibers, Phys. Rev. Lett. **126**, 023902 (2021).

[43] F. Mangini, M. Ferraro, M. Zitelli, V. Kalashnikov, A. Niang, T. Mansuryan, F. Frezza, A. Tonello, V. Couderc, A.B. Aceves, S. Wabnitz, Rainbow Archimedean spiral emission from optical fibres, Scientific Reports **11**, 13030 (2021).

[44] K. Stefanska, P. Béjot, K. Tarnowski, B. Kibler Experimental Observation of the Spontaneous Emission of a Space-Time Wavepacket in a Multimode Optical Fiber, ACS Photonics **10**, 727 (2023).

[45] L.G. Wright, S. Wabnitz, D. Christodoulides, F. Wise, Ultrabroadband dispersive radiation by spatiotemporal oscillation of multimode waves, Phys. Rev. Lett. **115**, 223902 (2015).

[46] K. Krupa, A. Tonello, A. Barthélémy, V. Couderc, B.M. Shalaby, A. Bendahmane, G. Millot, S. Wabnitz, Observation of geometric parametric instability induced by the periodic spatial self-imaging of multimode waves, Phys. Rev. Lett. **116**, 183901 (2016).

[47] K. Krupa, C. Louot, V. Couderc, M. Fabert, R. Guenard, B. M. Shalaby, A. Tonello, D. Pagnoux, P. Leproux, A. Bendahmane, R. Dupiol, G. Millot, and S. Wabnitz, Spatiotemporal characterization of supercontinuum extending from the visible to the mid-infrared in a multimode graded-index optical fiber, Opt. Lett. **41**, 5785 (2016).

[48] Z. Eslami, L. Salmela, A. Filipkowski, D. Pysz, M. Klimczak, R. Buczynski, J.M. Dudley, G. Genty, Two octave supercontinuum generation in a non-silica graded-index multimode fiber, Nature Comm. **13**, 2126 (2022).

[49] A. Chiocchetta, P.E. Larré, I. Carusotto, Thermalization and Bose-Einstein condensation of quantum light in bulk nonlinear media, Europhys. Lett. **115**, 24002 (2016).

[50] M.J. Davis, S.A. Morgan, K. Burnett, Simulations of Bose Fields at Finite Temperature, Phys. Rev. Lett. **87**, 160402 (2001).

[51] C. Connaughton, C. Josserand, A. Picozzi, Y. Pomeau, S. Rica, Condensation of classical nonlinear waves, Phys. Rev. Lett. **95**, 263901 (2005).

[52] P.B. Blakie, A.S. Bradley, M.J. Davis, R.J. Ballagh, C.W. Gardiner, Dynamics and statistical mechanics of ultra-cold Bose gases using c-field techniques, Advance in Physics **57**, 363-455 (2008).

[53] S. Nazarenko, *Wave Turbulence* (Springer, Lectures Notes in Physics, 2011).

[54] G. Krstulovic, M. Brachet, Energy cascade with small-scale thermalization, counterflow metastability, and anomalous velocity of vortex rings in Fourier truncated Gross-Pitaevskii equation, Phys. Rev. E **83**, 066311 (2011).

[55] N. Santic, A. Fusaro, S. Salem, J. Garnier, A. Picozzi, R. Kaiser, Nonequilibrium precondensation of classical waves in two dimensions propagating through atomic vapors, Phys. Rev. Lett. **120**, 055301 (2018).

[56] Y. Castin, R. Dum, E. Mandonnet, A. Minguzzi and I. Carusotto, Coherence properties of a continuous atom
laser, Journal of Modern Optics **47**, 2671 (2000).

[57] E. Mandonnet, A. Minguzzi, R. Dum, I. Carusotto, Y. Castin, J. Dalibard, Evaporative cooling of an atomic beam, Eur. Phys. J. D **10**, 9 (2000).

[58] T. Lahaye, Z. Wang, G. Reinaudi, S.P. Rath, J. Dalibard, D. Guéry-Odelin, Evaporative cooling of a guided rubidium atomic beam, Phys. Rev. A **72**, 033411 (2005).

[59] S.E. Olson, R.R. Mhaskar, G. Raithel, Continuous propagation and energy filtering of a cold atomic beam in a long high-gradient magnetic atom guide, Phys. Rev. A **73**, 033622 (2006).

[60] C.-C. Chen, S. Bennetts, R.G. Escudero, B. Pasquiou, F. Schreck, Continuous guided strontium beam with high phase-space density, Phys. Rev. Appl. **12**, 044014 (2019).

[61] V.S. L'vov, S.V. Nazarenko, Discrete and mesoscopic regimes of finite-size wave turbulence, Phys. Rev. E **82**, 056322 (2010).

[62] M. Onorato, Y.V. Lvov, G. Dematteis, S. Chibbaro, Wave Turbulence and thermalization in one-dimensional chains, Physics Reports **1040**, 1 (2023).

[63] M. Kolesik, J.V. Moloney, and M. Mlejnek, Unidirectional Optical Pulse Propagation Equation, Phys. Rev. Lett. **89**, 283902 (2002).

[64] A. Couairon and A. Mysyrowicz, Femtosecond filamentation in transparent media, Physics Reports **441**, 47 (2007).

[65] A. Couairon, E. Brambilla, T. Corti, D. Majus, O. de J. Ramirez-Gongora, and M. Kolesik, Practitioner's guide to laser pulse propagation models and simulation, Eur. Phys. J. Special Topics **199**, 5-76 (2011).

[66] A. D. Bandrauk, E. Lorin, and J.V. Moloney, Laser Filamentation, Mathematical Methods and Models (Springer, 2016).

[67] P. Béjot, Multimodal unidirectional pulse propagation equation, Phys. Rev. E **99**, 032217 (2019).

[68] V.E. Zakharov, V.S. L'vov, G. Falkovich, *Kolmogorov Spectra of Turbulence I* (Springer, Berlin, 1992).

[69] A.C. Newell, S. Nazarenko, L. Biven, Wave turbulence and intermittency, Physica D **152**, 520 (2001).

[70] A.C. Newell, B. Rumpf, Wave Turbulence, Annu. Rev. Fluid Mech. **43**, 59 (2011).

[71] S.K. Turitsyn, S.A. Babin, E.G. Turitsyna, G. E. Falkovich, E. Podivilov, D. Churkin, Optical Wave Turbulence, in *Advances in Wave Turbulence*, World Scientific Series on Nonlinear Science Series A, Vol. 83, edited by V.I. Shrira and S. Nazarenko (World Scientific, Singapore, 2013).

[72] D. Churkin, I. Kolokolov, E. Podivilov, I. Vatnik, S.

Vergeles, I. Terekhov, V. Lebedev, G. Falkovich, M. Nikulin, S. Babin, and S.K. Turitsyn, Wave kinetics of a random fibre laser, Nature Comm. **2**, 6214 (2015).

[73] J. Laurie, U. Bortolozzo, S. Nazarenko, S. Residori, One-dimensional optical wave turbulence: experiment and theory, Physics Reports **514**, 121-175 (2012).

[74] A. Picozzi, J. Garnier, T. Hansson, P. Suret, S. Randoux, G. Millot, D.N. Christodoulides, Optical wave turbulence: Toward a unified nonequilibrium thermodynamic formulation of statistical nonlinear optics, Physics Reports **542**, 1-132 (2014).

[75] Y. Zhu, B. Semisalov, G. Krstulovic, S. Nazarenko, Self-similar evolution of wave turbulence in Gross-Pitaevskii system, Phys. Rev. E **108**, 064207 (2023).

[76] Z. Wang, W. Fu, Y. Zhang, H. Zhao, Wave-turbulence origin of the instability of Anderson localization against many-body interactions, Phys. Rev. Lett. **124**, 186401 (2020).

[77] N. Cherroret, T. Scoquart, D. Delande, Coherent multiple scattering of out-of-equilibrium interacting Bose gases, Annals of Physics **435**, 168543 (2021).

[78] A.Y. Ramos, C. Shi, L.J. Fernandez-Alcazar, D.N. Christodoulides, T. Kottos, Theory of localization-hindered thermalization in nonlinear multimode photonics, Communications Physics **6**, 189 (2023).

[79] K.M. Frahm and D.L. Shepelyansky, Nonlinear Perturbation of Random Matrix Theory, Phys. Rev. Lett. **131**, 077201 (2023).

[80] P. Haldar, S. Mu, B. Georgeot, J. Gong, C. Miniatura, G. Lemarié, Rayleigh-Jeans prethermalization and wave condensation in a nonlinear disordered Floquet system, EPL **144**, 63001 (2024).

[81] S. Erne, R. Bücker, T. Gasenzer, J. Berges, and J. Schmiedmayer, Universal dynamics in an isolated one-dimensional Bose gas far from equilibrium, Nature **563**, 225 (2018).

[82] M. Kolesik, J. V. Moloney, Nonlinear optical pulse propagation simulation: From Maxwell's to unidirectional equations, Phys. Rev. E **70**, 036604 (2004).

[83] A. Biasi, O. Evnin, and B. A. Malomed, Fermi-Pasta-Ulam phenomena and persistent breathers in the harmonic trap, Phys. Rev. E **104**, 034210 (2021).

[84] E. Podivilov, D. Kharenko, V. Gonta, K. Krupa, O.S. Sidelnikov, S. Turitsyn, M.P. Fedoruk, S.A. Babin, S. Wabnitz, Hydrodynamic 2D turbulence and spatial beam condensation in multimode optical fibers, Phys. Rev. Lett. **122**, 103902 (2019).