# Wave turbulence, thermalization and multimode locking in optical fibers

M. Ferraro, K. Baudin, M. Gervaziev, A. Fusaro, Antonio Picozzi, J. Garnier, G. Millot, D. Kharenko, E. Podivilov, S. Babin, et al.

► To cite this version:

M. Ferraro, K. Baudin, M. Gervaziev, A. Fusaro, Antonio Picozzi, et al. Wave turbulence, thermalization and multimode locking in optical fibers. Physica D: Nonlinear Phenomena, 2025, 481, pp.134758. hal-05072552

HAL Id: hal-05072552

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

Submitted on 19 May 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
Wave turbulence, thermalization and multimode
locking in optical fibers

M. Ferraro<sup>a,b</sup>, K. Baudin<sup>c</sup>, M. Gervaziev<sup>d,e</sup>, A. Fusaro<sup>f</sup>, A. Picozzi<sup>c</sup>, J. Garnier<sup>g</sup>, G. Millot<sup>c</sup>, D. Kharenko<sup>d,e</sup>, E. Podivilov<sup>d,e</sup>, S. Babin<sup>d,e</sup>, F. Mangini<sup>h</sup>, S. Wabnitz<sup>b</sup>

<sup>a</sup>Department of Physics, University of Calabria, Via P. Bucci, Rende, 87036, (CS), Italy

<sup>b</sup>Department of Information Engineering, Electronics, and Telecommunications, Sapienza University of Rome, Via Eudossiana 18, Rome, 00184, Italy

<sup>c</sup>Université Bourgogne Europe, CNRS, Laboratoire Interdisciplinaire Carnot de Bourgogne ICB, UMR 6303, 21000 Dijon, France

<sup>d</sup>Department of Physics, Novosibirsk State University, Pirogova 1, Novosibirsk, 630090, Russia

<sup>e</sup>Institute of Automation and Electrometry, SB RAS, Novosibirsk, 630090, Russia

<sup>f</sup>Université Paris-Saclay, CEA, LMCE, Bruyères-le-Châtel, 91680, France

<sup>g</sup>CMAP, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, 91120, France

<sup>h</sup>Department of Engineering, Niccolò Cusano University, Rome, 00166, Itali

Abstract

We present a comprehensive overview of recent advances in theory and ex-
periments on complex light propagation phenomena in nonlinear multimode
fibers. On the basis of the wave turbulence theory, we derive kinetic equations
describing the out-of-equilibrium process of optical thermalization toward the
Rayleigh-Jeans (RJ) equilibrium distribution. Our theory is applied to explain
the effect of beam self-cleaning (BSC) in graded-index (GRIN) fibers, whereby
a speckled beam transforms into a bell-shaped beam at the fiber output as
the input peak power grows larger. Although the output beam is typically
dominated by the fundamental mode of the fiber, higher-order modes (HOMs)
cannot be fully depleted, as described by the turbulence cascades associated to
the conserved quantities. We theoretically explore the role of random refrac-
tive index fluctuations along the fiber core, and show how these imperfections
may turn out to assist the observation of BSC in a practical experimental set-
ting. This conclusion is supported by the derivation of wave turbulence kinetic
equations that account for the presence of a time-dependent disorder (random
mode coupling). The kinetic theory reveals that a weak disorder accelerates
the rate of RJ thermalization and beam cleaning condensation. On the other
hand, although strong disorder is expected to suppress wave condensation,
the kinetic equation reveals that an out-of-equilibrium process of condensa-
tion and RJ thermalization can occur in a regime where disorder predominates
over nonlinearity. In general, the kinetic equations are validated by numerical
simulations of the generalized nonlinear Schrodinger equation. We outline a
series of recent experiments, which permit to confirm the statistical mechanics
approach for describing beam propagation and thermalization. For example,
we highlight the demonstration of entropy growth, and point out that there are
inherent limits to peak-power scaling in multimode fiber lasers. We conclude by pointing out the experimental observation that BSC is accompanied by an effect of modal phase-locking. From the one hand this explains the observed preservation of the spatial coherence of the beam, but also it points to the need of extending current descriptions in future research.

Keywords:

Nonlinear optics, Statistical mechanics, Wave turbulence, Multimode fibers, Beam self-cleaning, Wave thermalization, Wave condensation

# 1 Introduction and theoretical framework

Up to now, the vast majority of optical fiber technologies are based on the use of singlemode fibers (SMFs) [1]. Historically, multimode fibers (MMFs) have been manufactured and studied well before SMFs: an interesting example is the first supercontinuum generation observation in an MMF, performed back in 1978 by by Lin et al. [2]. The reason for SMFs to outplace MMFs can be easily explained. Owing to their reduced core diameter, SMFs only support one mode, so that the spatial of the guided beam are fixed by the fiber geometry. Whereas an optical beam that is carried by an MMF is typically transported by hundreds or even thousands of spatial modes, whose relative occupancy may vary upon propagation owing to random linear mode coupling. Since each mode propagates with its own phase and group velocity, multimode interference at the output of a MMF will lead to highly unstable speckled intensity patterns with extremely low beam quality, which prevents their use for imaging applications that involve focused beams. Most importantly, for lightwave communications the use of MMFs is limited to very short distances only. This is because modal dispersion largely overcomes chromatic dispersion on the one hand, and also because its post-compensation by digital means (e.g., by multiple-input, multiple-output (MIMO) techniques) is very complex unless just a few modes are present.

Nevertheless, a renewed interest in light propagation in nonlinear MMFs has emerged over the past two decades [3]. This is mostly driven by the exponential growth in optical network traffic, and the potential of MMFs for increasing transmission capacity via spatial division (or mode division) multiplexing. In addition, the possibility offered by MMFs (thanks to their large cross-section, hence increased power threshold for optical damage) of permitting high-energy beam delivery is of utmost interest for scaling up the power generated by fiber lasers, and for multiphoton imaging applications.

In this respect, studying nonlinear beam propagation in MMF is very relevant [4]. Indeed, nonlinearity may be detrimental for both telecommunication applications and for high-power fiber laser sources, and many theoretical and experimental studies have been carried out in order to mitigate, and possibly suppress, the impact of nonlinear effects in MMFs. In contrast, nonlinear effects may be beneficial for the development of supercontinuum light sources,
as well as for achieving high-resolution biomedical imaging, e.g., in the field of endoscopy.

Still, the fundamental issue of how to deal with being able of theoretically predicting the outcome of the complex nonlinear interaction among such a large number of modes, i.e., of degrees of freedom, holds. As a matter of fact, in the presence of third-order nonlinear effects, such as the Kerr effect, Raman and Brillouin scattering, the control of mode interactions can be achieved in ways that have no counterpart in the singlemode fiber case [5, 6]. Some examples of robust wave structures which emerge from complex nonlinear mode mixing in MMFs are multimode optical solitons [7] and ultra-wideband parametric sidebands [8] [9, 10].

Looking backwards in time, a similar complex problem was faced by the physicists of the 19th century, when addressing the description of the physical properties of a gas in the framework of purely classical mechanics. Indeed, an ideal gas is nothing more than a system of N particles, whose motion solves Newton's equation. As it is well known, a dynamical system with three interacting bodies may be chaotic. Of course, whenever $N \gg 1$, solving the problem of motion by brute force numerical integration of a system of N differential equations may become virtually impossible. In a gas of classical particles, N is of the order of the Avogadro's number ($\sim 10^{23}$). With no supercomputer at hand, Boltzmann and his fellow 1800s physicists had to find a clever way to conciliate classical mechanics with thermodynamics, without spending their lifetime to solve a complex system of coupled equations. Such a clever way was the development of statistical mechanics: the macroscopic properties of a mechanical system with many degrees of freedom are obtained by averaging over a statistical ensemble of the microscopic properties.

In this way, one does not need to know the exact position and velocity of all particles in a gas, for describing its macroscopic state. Only average (i.e., thermodynamic) quantities matter. These are the temperature, the pressure, the chemical potential, the entropy, and so on. A similar approach can be used for nonlinear multimode photonics.

## 1.1 The statistical mechanics theory

One of main successes of the statistical mechanics theory was to provide a kinetic theory of gases. According to this theory, at thermal equilibrium (which is defined as the state of maximum entropy) the velocities of all particles, which are supposed to be distinguishable, follow a Maxwell-Boltzmann probability distribution. Applying a statistical mechanics approach to nonlinear multimode *classical* waves leads to a different statistics for the mode power distribution. According to the theory, as long as the nonlinearity can be considered *weak* and the number of modes is finite [11], the particles that compose a guided optical beam (i.e., the photons) follow a Rayleigh-Jeans (RJ) law. This can be easily understood by the following argument. Photons are bosons, which are indistinguishable and obey the Bose-Einstein distribution. This boils down to the RJ distribution in the classical limit, i.e., whenever the number of photons per mode is sufficiently high. In this review, we disregard the quantum
nature of light, and only consider the classical limit in which the modes of the
waveguide are highly populated. As such, we will always refer to the RJ law as
the thermal equilibrium distribution.

Interestingly, the statistical mechanics approach describing nonlinear beam
propagation in MMF has been developed by following two different approaches.
First, in 2011, a weak wave-turbulence theory was introduced by Ascheri et al..
Their work involved the processes of classical wave thermalization and conden-
sation in GRIN MMFs [11], whose experimental verification was only carried
out nine years later by Baudin et al. [12]. On the other hand, Wu et al. and Parto
et al., proposed in 2019 a thermodynamic theory of highly multimode nonlin-
ear optical systems, developing an equation of state at thermal equilibrium that
links together all of the thermodynamic parameters [13, 14]. The results of that
work, mostly based on probabilistic (i.e., entropic) concepts, were then refined
by Makris et al., who developed a theoretical framework for describing the
thermalization process as a consequence of photon-photon interactions [15], in
analogy with the textbook statistical mechanics of ideal gases. The statistical
mechanics approach has been subsequently applied and extended to various
different forms of multimode optical systems [16, 17, 18], e.g., to study the
thermalization in non-Hermitian optical lattices [19], or to refine the notion of
optical pressure [20, 17, 21, 22, 23, 24].

As long as equilibrium states are involved, the two approaches provide
exactly the same results, i.e., they lead to the same RJ distribution. This is
because the two theories rely on the very same pillars. Both theories consider
a weakly nonlinear fiber supporting a finite number of guided modes; upon
propagation, an optical beam conserves both its power (number of photons)
and "momentum", i.e., the system's Hamiltonian (referred to either as energy
[11, 25] or as internal energy [13, 15]). As a matter of fact, in standard cylindrical
MMF there is a third conserved quantity, i.e., the orbital angular momentum
of light [26, 27]. However, for the sake of simplicity, in this review we will
limit ourselves to consider optical beams that carry no orbital angular momen-
tum. The beam energy has two contributions, which can be labeled either
linear and nonlinear, or kinetic and potential. In any case, in the weakly non-
linear regime, i.e., whenever the power of the optical beam is sufficiently low,
the potential energy (nonlinearity) can be neglected in comparison with the
kinetic energy. Each conserved quantity is associated with a thermodynamic
parameter. Specifically, the conservation of the number of photons and of the
energy lead to the definition of an optical chemical potential and an optical
temperature, respectively. Notice that statistical mechanics has been also ex-
tended beyond the weakly nonlinear regime, by computing elaborate partition
functions [28, 29, 30, 31, 32, 33]. Along this way, the statistical properties of in-
coherent optical fields have been also analyzed by applying methods inherited
from complex systems and spin-glass theory [34, 35, 36, 37].

The true potency of the statistical mechanics approach for predicting the
evolution of weakly nonlinear multimode waves is its simplicity and effective-
ness: by a trivial algebraic calculation, one obtains the (average) spatial profile
of a guided beam into an MMF, based on the mere knowledge of its injection
conditions. This is the main reason behind the outburst of theoretical and exper-
imental studies on this topic which appeared in recent literature. Without the
statistical mechanics approach, determining the beam at the output of a highly
multimode fiber would require solving the nonlinear wave equation numeri-
cally. With standard numerical tools, this may last even days of computation.
The trick behind the statistical mechanics approach is somehow equivalent to
"measuring" the temperature and chemical potential of the beam for determin-
ing its properties. Analogously to Boltzmann's approach to classical particles
gases, one does not need to know the amplitude and phase of all individual
modes in order to know, for instance, its temperature. It has to be mentioned,
however, that, at variance with classical thermodynamics, there is no ther-
mometer available here to measure the beam temperature when dealing with
optical beams. The optical temperature has a pure statistical meaning, which
can even take negative values, owing to the finite energy spectrum of multi-
mode optical systems, as originally predicted in Refs.[13, 14] and subsequently
experimentally reported [38, 39]. In a loose sense, the optical temperature re-
flects the amount of randomness of a speckled beam at equilibrium, and it can
be computed by the sole knowledge of the beam's conserved quantities (kinetic
energy and number of particles), which univocally determine the equilibrium
mode power distribution.

At this point, one may object that the thermodynamic equivalence picture
that we have discussed only strictly holds at thermal equilibrium. Therefore,
one may naturally wonder whether there is a way, aside from equilibrium statis-
tical mechanics approaches, to provide a detailed nonequilibrium description
of the process of thermalization in the weakly nonlinear regime. This is indeed
possible by recurring to the help of wave turbulence theory.

1.1.1 *The wave turbulence theory*

The wave turbulence theory occupies a rather special place on the roadmap
of modern science, at the interface between applied mathematics, fluid dy-
namics, statistical physics, and engineering [40, 41, 42, 43]. It has potential
applications in a diverse range of subjects, including oceanography, plasma
physics, and condensed matter physics. Optical wave turbulence also consti-
tutes a growing field of research covering various topics in modern optics, e.g.,
superfluid turbulence [44, 45], rogue wave generation [46, 47, 48, 49], super-
continuum generation [50, 51, 52, 53, 54], integrable turbulence [55, 56], light
condensation [57, 11, 58, 59, 60, 61] and thermalization [62, 63, 64], nonequi-
librium transport in multimode optical systems [65, 66], or dissipative cavities
and lasers [67, 68, 69, 70, 71] – see Refs. [72, 73] for reviews.

Originally, the wave turbulence theory has been the subject of intense inves-
tigations in the context of plasma physics [74], in which it is often referred to the
so-called "random phase-approximation" approach [40, 74, 75, 76, 77, 78]. This
approach may be considered as a convenient way of interpreting the results of
the more rigorous technique based on a multi-scale expansion of the cumulants
of the nonlinear field, as originally formulated in Refs. [79, 80, 81]. This theory
has been reviewed in Ref. [82], and studied in more detail through the analysis
of the probability distribution function of the random field in Refs.[41]. In a
loose sense, the so-called 'random phase approximation' may be considered
as being justified when phase information becomes irrelevant to the wave in-
teraction, due to the strong tendency of the waves to decohere. The random
phases can thus be averaged out, to obtain a weak turbulence description of
the incoherent wave interaction, which is formally based on irreversible kinetic
equations [40]. It turns out that, in spite of the formal reversibility of the equa-
tion governing wave propagation, the kinetic equation describes an irreversible
evolution of the field towards its thermodynamic equilibrium. Such an equilib-
rium state refers to the Rayleigh-Jeans spectrum, whose tails are characterized
by an equipartition of energy among the modes. The mathematical statement of
such an irreversible process is the *H*-theorem of entropy growth, whose origin
is analogous to the Boltzmann's *H*-theorem relevant for gas kinetics [83].

Referring back to light propagation in GRIN MMFs, the wave turbulence
theory is capable of describing the mode distribution associated with the so-
called beam self-cleaning (BSC) effect, a striking phenomenon that has been
attracting a great deal of interest over the last few years.

1.2 *The beam self-cleaning effect*

First reported by Krupa et al. in 2016 [84, 85] for nanosecond pulses and
later confirmed for femtosecond pulses by Lui et al. [86], the BSC effect in GRIN
MMFs refers to the phenomenon whereby an input beam evolves, during its
propagation through complex and spatially disordered states, to form a stable,
well-defined output beam with a typical near-Gaussian profile. Such an effect
can be numerically reproduced by using different approaches: e.g., the full-
field model [85], as well as the generalized multimode nonlinear Schrödinger
equations (GMM-NLSE) [87]. Within the thermodynamic framework, the oc-
currence of BSC is seen as the establishment of a state of thermal equilibrium
[88], as depicted in Fig. 1. A thermalized output beam, usually referred to as
"clean", turns out to be remarkably stable versus external perturbations. This
has paved the way for the possible development of a novel class of MMF-based
applications. As previously discussed, most of current fiber technologies are
based on SMFs. Since SMFs only support one mode, they inherently preserve
the beam's spatial quality upon its propagation, even under the action of ex-
ternal perturbations, such as fiber banding and squeezing. When it comes to
MMFs, the multitude of fiber modes that can interfere with each other do not
ensure the preservation of the beam spatial profile, even in the absence of any
onlinear effects. Moreover, the output beam may significantly vary if external
perturbations are applied, e.g. temperature fluctuations or mechanical stress.

Indeed, preserving the spatial quality of the beam is a key feature when
transporting it via optical fibers. For instance, who would buy a laser whose
output beam is speckled and variable in time because of weak environmental
variations? As a result of this rhetorical question, until a few years back most of
the research on optical fibers was focused on SMFs. It is in this framework that
the discovery of BSC was groundbreaking. Guaranteeing a high beam quality
together with high-energy pulse delivery allows for scaling up the power of SMF
lasers, as demonstrated by Teğin et al. [89] and for enhancing the resolution of
multiphoton imaging devices as shown by Moussa et al. [90, 91].

<figure><img src="image_2.png" /><figcaption>Figure 1: Toy figure of BSC as a result of classical wave thermalization/condensation.</figcaption></figure>

1.3 *The goal and the content of this review*

In essence, we aim to offer a comprehensive overview of nonlinear multi-mode fiber optics within a statistical mechanics framework. In particular, we focus on the wave turbulence approach, since this allows for capturing the whole process of thermalization, i.e., the evolution of a multimode laser beam through intermediate non-equilibrium states, eventually reaching RJ equilibrium. The wave turbulence theory specifically allows for the characterization of beam self-cleaning by providing a characteristic propagation length scale needed to achieve thermal equilibrium. This is important, in particular, when one considers the key role played in the thermalization process by random fluctuations of the MMFs' refractive index profile due to fabrication imperfections.

We give particular emphasis to the case of BSC in GRIN MMFs. From the fundamental point of view, since its first observations it was clear that observing BSC requires the combined action of the Kerr and the spatial self-imaging (SI) effects [92]. SI arises from the unique parabolic refractive index profile of GRIN MMFs (BSC has not yet been clearly observed to occur in step-index fibers). These cause the modes to mix and redistribute their energy, eventually favoring an energy transfer towards the fundamental mode. The outcome of BSC is that the fundamental mode soaks up as much power as possible from the high-order modes (HOMs). However, according to the statistical mechanics theory, HOMs cannot be completely emptied, because this would result in a breach of the conservation laws of light power and kinetic energy. The fact that a BSC beam is not composed by the sole fundamental mode was experimentally verified by Leventoux in 2020 [93]. This excluded possible explanations of BSC as a result of mode-dependent losses, thus strengthening the validity of the statistical mechanics picture, which, in contrast, is conservative. Indeed, the
conservation laws drive the nonlinear beam dynamics in a way that thermal
equilibrium is established, as a result of optical entropy maximization.

In this work, we largely focus on the role played by random disorder, which
has been shown to facilitate the BSC process [94]. We frame the role played
by random disorder within wave turbulence theory. As a matter of fact, the
wave turbulence theory has been essentially developed by ignoring the impact
of disordered fluctuations of the propagating nonlinear medium [40, 42, 41].
However, the evolution of a turbulent wave system can be heavily affected by
environmental disorder. This plays a fundamental role in ocean waves, where
the scattering by disorder can originate from random variations of the bottom
topography in oceans [95]. In fiber optics, disorder originates from random
fluctuations of the refractive index of the core material. In particular, disorder
is known to impact light propagation in MMFs, due to imperfections and en-
vironmental perturbations, a feature that is relevant for instance to endoscopic
imaging [96]. So far, the impact of disorder on turbulent flows has been essen-
tially studied by means of numerical simulations, while only few theoretical
developments have considered purely time-independent disordered fluctua-
tions [97, 98, 99, 100, 101, 102]. On the other hand, light propagation in MMFs
is characterized by random refractive index fluctuations that vary along the
propagation, and are thus ‘time’-dependent. In Sec. 2 we review recent works
[103, 25, 104], in which the wave turbulence theory has been extended to ac-
count for the presence of a time-dependent disorder, which is inherent to light
propagation in MMFs.

In Sec. 3, we present a selection of experimental investigations aimed at
confirming the predictions of the statistical mechanics approach for nonlinear
beam propagation in GRIN MMFs. We first consider the notable case of BSC.
The validity of the statistical mechanics approach to BSC was first proved by the
experiments of Baudin et al. in 2020 [12]. Then, thanks to the development of
a field-projection mode decomposition (MD) based of self-referenced interfer-
ometric technique [105] and holographic MD [106] tools, a direct measurement
of the RJ mode distribution was reported by Pourbeyram et al. [107] and, later,
by Mangini et al. [108] in 2022. For a review of these early works, we refer the
reader to Ref.[88]. Here, instead, we focus on more recent advances. In par-
ticular, we highlight experimental evidences of the entropy growth during the
thermalization of either one [109] or two waves [110]. Next we discuss recent
intriguing experiments: they demonstrate that, by increasing the temperature
of a self-cleaned (i.e., at equilibrium) beam, one observe a progressive loss its
beam quality. This result establishes a practical limit to the power up-scaling
of MMF lasers. Then, we focus on an aspect which has been largely debated
within the optical statistical mechanics community, i.e., the framing of the sta-
tistical ensemble in experiments. Indeed, one of the pillars of the statistical
mechanics theory is the averaging of physical observables over a meaningful
statistical ensemble. As aforementioned, in the case of optical thermodynamics
or wave turbulence, one generally averages over the mode phases. However,
virtually all the BSC experimental observations did not involve a properly said
phase average. In this sense, the only exception reported in the literature is
Ref. [12], where an average over an ensemble of injected speckle beams was considered. At last, we present the recent observation of a modal phase-locking mechanism that accompanies BSC [111].

Finally, in Sec. 4, we provide our view on the BSC effect in light of recent results that have questioned its interpretation as a pure thermalization process. Then we conclude by tracing what will be, in our opinion, the foreseeable progress of research in the domain of the statistical mechanics of nonlinear MMF optics.

## 2 Wave turbulence in multimode fibers with disorder

### 2.1 Basic principles and conservation laws

Before introducing the technical formalism inherent to the wave turbulence theory, we aim to provide some simple physical insight into light thermalization in MMFs, grounded on basic principles and conservation laws. When considering the purely 2D spatial evolution of the optical field in a MMF, we may expand the monochromatic field amplitude $\psi(z, \mathbf{r}) = \sum_{\ell,m} a_{\ell,m}(z) u_{\ell,m}(\mathbf{r}) \exp(i\beta_{\ell,m}z)$ on the orthonormal basis of the eigenmodes $u_{\ell,m}(\mathbf{r}) = R_{\ell,m}(r) \exp(im\phi)$ of the waveguide, where $R_{\ell,m}(r)$ are the radial components of the eigenmodes, with $z$ the spatial coordinate along the fiber axis and $\mathbf{r} = (r, \phi)$ the transverse coordinate. Here, $\beta_{\ell,m}$ denote the propagation constants [with $\beta_{\ell,m} = \beta_0(2\ell + |m|)$ for the case of the parabolic graded-index fiber], while $|a_{\ell,m}|^2$ denote the power occupancy of the optical field, $(\ell, m)$ being the radial and azimuthal numbers (with $\ell = 0, 1, 2, ...$ and $m = 0, \pm 1, \pm 2, ...$).

The optical field is characterized by the conservation of different quantities during its propagation through the fiber. The optical power, $N = \sum_{\ell,m} |a_{\ell,m}|^2$ (or 'wave-action' in the wave turbulence terminology, or particle number in a corpuscular picture) is a conserved quantity, as well as the orbital angular momentum, $M = \sum_{\ell,m} m |a_{\ell,m}|^2$. Furthermore, experiments demonstrating RJ thermalization in MMFs have thus far been conducted in a weakly nonlinear regime. Accordingly, the (kinetic) energy $E = \sum_{\ell,m} \beta_{\ell,m} |a_{\ell,m}|^2$, that is the linear contribution to the Hamiltonian, is also conserved through the propagation.

All conserved quantities ($N$, $E$, $M$) are completely determined by the injection conditions of the laser beam into the MMF [112]. During the propagation in the MMF, Kerr nonlinearity leads to energy exchanges among the fiber modes, similarly to particle collisions in a classical gas, thus shuffling the occupation of the fiber modes, $a_{\ell,m}(z)$. One could expect that the optical speckled beam reaches thermodynamic equilibrium over a finite "time" (i.e., a finite propagation length $z$), owing to the well-known Onsager's principle of detailed balance [113]. When applied to mode interactions, this principle implies that at thermal equilibrium (the thermalized state of an optical beam), every elementary four-wave mixing process that transfers energy among the modes is equally probable as its reverse process. At equilibrium, the average power occupation in a specific mode, $n_{\ell,m}^{eq} = \langle |a_{\ell,m}|^2 \rangle$, follows a RJ distribution. The
detailed expression of the equilibrium distribution is obtained by maximizing the optical entropy, whose appropriate expression, $S = \sum_{\ell,m} \log(n_{\ell,m})$, can be derived from a statistical mechanics approach [13, 15], or alternatively from a H-theorem of entropy growth [11], as discussed below in the framework of the wave turbulence theory, see section 2.3. Maximizing the entropy $S[n_{\ell,m}]$ subject to the constraints of the conservation of the power $N = \sum_{\ell,m} n_{\ell,m}$, the energy $E = \sum_{\ell,m} \beta_{\ell,m} n_{\ell,m}$, and the angular momentum $M = \sum_{\ell,m} mn_{\ell,m}$, readily gives [112, 27]:

$$n_{\ell,m}^{eq} = \frac{1}{a + b\beta_{\ell,m} + cm}, \qquad (1)$$

where $a, b$ and $c$ are three Lagrange's multipliers associated to the three conserved quantities, $N, E$ and $M$.

The conservation of angular momentum provided unforeseen insights into the physical mechanisms underlying spatial beam cleaning. For instance, an inspection of Eq.(1) reveals that the fundamental mode of the fiber is not necessarily the most populated in the presence of a non-vanishing angular momentum [27], a remarkable feature that has been demonstrated experimentally in Ref.[112].

It should be noted that launching a beam carrying a non-zero angular momentum in a MMF requires specific injection conditions. In this section, we consider the general framework of the wave turbulence approach, and typically assume initial conditions characterized by a random phase among the modes. More precisely, considering an average over a statistical ensemble of speckled beams with random phases, the average angular momentum vanishes at the fiber input, $M = \sum_{\ell,m} m \langle |a_{\ell,m}|^2 \rangle = 0$, which implies $c = 0$ in the RJ equilibrium distribution (1). In this way, one usually poses $a = -\mu/T$ and $b = 1/T$ following the standard procedure, where $\mu$ and $T$ are the chemical potential and the temperature in analogy to thermodynamics. The equilibrium RJ distribution then recovers the standard form:

$$n_{\ell,m}^{eq} = \frac{T}{\beta_{\ell,m} - \mu}. \qquad (2)$$

The RJ distribution determines the equilibrium properties of the optical field. However, the equilibrium approach discussed here provides no information about the characteristic "time" scales (i.e., propagation lengths) needed to reach such an equilibrium state. In particular, random mode coupling can lead to power and energy exchanges among the modes, which strongly affect the rate of thermalization to equilibrium, as described by the wave turbulence kinetic theory.

## 2.2 Impact of random mode coupling

As discussed in the introduction section, the wave turbulence theory has been developed by ignoring the presence of a structural disorder in the system [40, 82, 42, 78, 41, 72, 73]. So far, the impact of disorder on turbulent flows has been essentially studied by means of numerical simulations, while only
few theoretical developments have considered purely time-independent (i.e., $z$-independent) disordered fluctuations [97, 98, 99, 100, 101, 102]. On the other hand, light propagation in MMFs is known to be affected by random fluctuations of the longitudinal and transverse profiles of the index of refraction, as a consequence of external factors such as bending, twisting, tensions, kinks, or core-size variations in the fabrication process of the fiber. Such multiple physical perturbations introduce random polarization fluctuations as well as random coupling among the modes of the fiber. The specific mechanisms and models that describe how fiber imperfections impact light propagation in MMFs still remain an active topic of research [114, 115, 116, 117, 118, 119].

One usually assumes that for relative short propagation lengths (few dozen meters), the contribution of noise essentially arises from 'weak' disorder, which originates from polarization random coupling and random coupling among degenerate fiber modes [118, 120]. The presence of polarization random coupling does not prevent the conservation of the three conserved quantities ($N, E, M$), so that the expected RJ equilibrium has the form given in Eq.(1), while random coupling among degenerate modes modifies the angular momentum, so that the expected RJ distribution has the form (2). On the other hand, the contribution of 'strong' disorder arises for larger propagation lengths (few hundred meters) from random coupling among non-degenerate modes [118, 120]. The presence of strong disorder prevents the conservation of energy and angular momentum, so that the expected RJ equilibrium should reduce to an equipartition of the power among the modes, i.e., $b = c = 0$ in Eq.(1), resulting in a parabolic intensity profile in graded-index fibers, and a rectangular intensity profile in step-index fibers [121, 122]. In the following sections, we derive the kinetic equation for the case of weak disorder by considering polarization random fluctuations in Sec. 2.3, while the case of strong disorder will be addressed in Sec. 2.4. Unexpectedly, the kinetic theory reveals that thermalization to the RJ equilibrium (2) may even occur in the presence of strong disorder, see Sec. 2.4.

The impact of weak disorder in graded-index MMFs is also known to play an important role in the process of thermalization [94, 103, 25]. In particular, it was shown in Ref.[25] that in the absence of disorder, a strong correlation among the modes is preserved for large propagation lengths, which leads to a phase-sensitive "coherent regime" of modal interaction. In this coherent regime, the modes experience a quasi-reversible exchange of power with each other, which tends to freeze the process of RJ thermalization, as illustrated in Fig. 2. Note that such an oscillatory behavior of the modal components, and the slowing down of RJ thermalization, are related to the existence of Fermi-Pasta-Ulam recurrences, as recently discussed in Ref.[123, 124] in the framework of the weakly nonlinear regime of the 2D NLS equation with a parabolic trapping potential. We will see that the source of noise introduced by weak disorder (polarization random fluctuations) is sufficient to break the coherent regime of modal interaction. The resulting random phase dynamics then leads to a turbulent incoherent regime of the modal components, which can be described in the framework of the wave turbulence theory.
## 2.3 Kinetic equation with weak disorder

Let us consider the case of weak disorder in the framework of polarization random fluctuations, and discuss its impact on the wave turbulence kinetic equation in MMFs. We first treat the conventional "continuous" wave turbulence regime that may be relevant for very highly multimode fibers, and then the "discrete" wave turbulence regime, relevant to typical experiments of beam cleaning. On the basis of recent developments on finite size effects in wave turbulence [125, 126, 127, 128, 129], a detailed derivation of the discrete wave turbulence kinetic equation accounting for weak disorder is provided in Section 2.3.5. We first consider a model of weak disorder, in which the modes experience an independent (decorrelated) noise. We note that, in the opposite limit where the modes experience the same disorder, the analysis reveals that the optical beam recovers a quasi-coherent regime of mode interactions, similar to that obtained in the absence of disorder. In other words, the fully mode-correlated model of disorder appears ineffective, and in this limit the optical beam does not exhibit a fast process of condensation. In this way, an intermediate model of partially correlated disorder is considered, where only degenerate modes experience the same noise. The theory reveals that this partially mode-correlated disorder is sufficient to re-establish an efficient acceleration of condensation.

### 2.3.1 Model equation

Let us consider the weak disorder regime, taking into account the main contribution of polarization random fluctuations [116, 117, 120]. Accordingly, we resort to a vector model of light propagation in a GRIN MMF that includes polarization mode dynamics. In the rectilinear polarization basis, the evolutions of the modal components $a_p(z) = (a_{p,x}, a_{p,y})^T$ are governed by

$$i \partial_z a_p = \beta_p a_p + \mathbf{D}_p(z) a_p - \gamma \mathbf{P}_p(\mathbf{a}), \quad (3)$$

with the nonlinear term

$$\mathbf{P}_p(\mathbf{a}) = \sum_{l,m,n} S_{plmn} \left( \frac{1}{3} \mathbf{a}_l^T \mathbf{a}_m \mathbf{a}_n^* + \frac{2}{3} \mathbf{a}_n^\dagger \mathbf{a}_m \mathbf{a}_l \right), \quad (4)$$

where $\beta_p$ are the eigenvalues verifying $\beta_p u_p(\mathbf{r}) = -\alpha \nabla^2 u_p(\mathbf{r}) + V(\mathbf{r}) u_p(\mathbf{r})$, with $\alpha = 1/(2n_{co}k_0)$, and $k_0 = 2\pi/\lambda$, $\lambda$ being the laser wavelength. We assume that the refractive index potential is parabolic-shaped, $V(\mathbf{r}) = q|\mathbf{r}|^2$ for $|\mathbf{r}| \le R$, $R$ being the radius of fiber core, and $q = k_0(n_{co}^2 - n_{cl}^2)/(2n_{co}R^2)$ with $n_{co}$ and $n_{cl}$ the refractive indexes of the core and the cladding, respectively. Accordingly, the eigenvalues of the parabolic potential read $\beta_p = \beta_0(p_x + p_y + 1)$, where $p$ labels the two integers $(p_x, p_y)$ that specify a mode in the Hermite-Gauss basis, and $\beta_0 = 2\sqrt{\alpha q}$. In Eq.(4), $S_{lmnp} = \int u_l(\mathbf{r}) u_m(\mathbf{r}) u_n^*(\mathbf{r}) u_p^*(\mathbf{r}) d\mathbf{r} / \int |u_0|^4(\mathbf{r}) d\mathbf{r}$, so that $S_{0000} = 1$, and $\gamma = k_0 n_2 / A_{eff}^0$, where $A_{eff}^0 = 1 / \int |u_0|^4(\mathbf{r}) d\mathbf{r}$ is the effective area of the fundamental mode, and $n_2$ the nonlinear Kerr coefficient. The sum in (4) is carried over the mode indices $l = (l_x, l_y)$, $m = (m_x, m_y)$, $n = (n_x, n_y)$
<figure><img src="image_3.png" /><figcaption>Figure 2: Coherent regime of mode interaction without disorder: Numerical simulations of the NLS Eq.(3) showing the evolution of the modal components $n_p(z) = |a_p(z)|^2$ (in the absence of disorder $D_p = 0$), starting from a coherent initial condition (a), starting from a 'speckle' beam with random phase among the modes (b). Evolutions of the modal populations $n_p/N$, fundamental mode $p=0$ (red dashed), $p=1$ (dark blue solid), $p=2$ (blue solid), $p=3$ (light blue solid), $p=4$ (cyan solid), $p=5$ (light green solid), $p=6$ (green solid), $p=7$ (yellow solid), $p=8$ (orange solid). The system exhibits a persistent oscillatory dynamics among the low-order modes. At complete thermal equilibrium the system would reach the condensate fraction $n_0^{eq}/N \simeq 0.68$ in (a), and $n_0^{eq}/N \simeq 0.58$ in (b). The power is $N = 47.5$ kW, $M = 120$ modes, $R = 26$ µm. Source: From Ref.[25].</figcaption></figure>
of the M propagative (or guided) modes. The superscripts *, T, † stand for the conjugate, the transpose, and the conjugate transpose operations.

The characteristic nonlinear length is $L_{nl} \sim 1/(\gamma N)$, while the linear characteristic length from Eq.(3) is $L_{lin} = 1/\beta_0$. In usual experiments of beam self-cleaning in GRIN fibers we can have $\beta_0 \sim 5 \times 10^3 \text{ m}^{-1}$, so that the optical beam evolves in the weakly nonlinear regime

$$L_{lin} \ll L_{nl}, \qquad (5)$$

where $L_{lin} \lesssim 1$ mm, while $L_{nl}$ is typically larger than 10 cm. As will be discussed below in more details, the weakly nonlinear regime of the experiments allows us to derive a kinetic equation in the framework of the weak turbulence theory.

Let us recall that longitudinal and transverse fluctuations in the refractive index profile of the MMF lead to a random coupling among the modes, as expressed by the coupled-mode theory [119, 116]. However, by following this procedure one obtains a scalar perturbation that does not introduce coupling among the polarization components, a feature commented in particular in Ref.[117]. For this reason, random polarization fluctuations are usually introduced in a phenomenological way, see e.g.[116, 117, 120]. Here, we consider the most general form of polarization disorder that conserves the power $N = \sum_p |a_p|^2$, and the (kinetic) energy $E = \sum_p \beta_p |a_p|^2$, of the optical beam. For this purpose, the Hermitian matrices $D_p(z)$ are expanded into the Pauli matrices, which are known to form a basis for the vector space of $2 \times 2$ Hermitian matrices. The matrices then have the form

$$D_p(z) = \sum_{j=0}^{3} \nu_{p,j}(z) \sigma_j, \qquad (6)$$

where $\sigma_j$ ($j = 1, 2, 3$) are the Pauli matrices and $\sigma_0$ is the identity matrix. The functions $\nu_{p,j}(z)$ are independent and identically distributed Gaussian real-valued random processes, with

$$\langle \nu_{p,j}(z) \nu_{p',j'}(z') \rangle = \sigma_{\beta}^2 \delta_{pp'}^{K} \delta_{jj'}^{K} \mathcal{R} \left( \frac{z-z'}{l_{\beta}} \right). \qquad (7)$$

Here $l_\beta$ is the correlation length of the random process and $\sigma_\beta^2$ is its variance. The normalized correlation function $\mathcal{R}$ is such that $\mathcal{R}(0) = O(1)$, $\int_{-\infty}^{\infty} \mathcal{R}(\zeta) d\zeta = O(1)$, and $\mathcal{R}(z/l_\beta) \to l_\beta \delta(z)$ in the limit $l_\beta \to 0$. We will consider Ornstein-Uhlenbeck processes:

$$d\nu_{p,j} = -\frac{1}{l_{\beta}} \nu_{p,j} dz + \frac{\sigma_{\beta}}{\sqrt{l_{\beta}}} dW_{p,j}(z)$$

where $W_{p,j}$ are independent Brownian motions. This means that $\nu_{p,j}$ are Gaussian processes with mean zero and covariance function of the form (7)
with $\mathcal{R}(\zeta) = \exp(-|\zeta|)/2$. We introduce an effective parameter of disorder $\Delta\beta = \sigma_{\beta}^2 l_{\beta}$ and the associated length scale

$$L_d = 1/\Delta\beta. \qquad (8)$$

We assume that disorder effects dominate nonlinear effects

$$L_d \ll L_{nl}, \qquad (9)$$

and that $l_{\beta} \ll L_d$ (or equivalently $\sigma_{\beta}l_{\beta} \ll 1$). Let us recall that usual experiments of beam self-cleaning are carried out in the weakly nonlinear regime [see Eq.(5)], so that $\beta_0^{-1} \ll L_d$ (or $\Delta\beta \ll \beta_0$). Aside from beam-cleaning experiments, in the following we will also consider a regime where the mode-spacing can be very small $\beta_0 \ll \Delta\beta$ (MMFs with huge number of modes $M \gg 1$), in order to address a regime described by a continuous wave turbulence approach. Finally, note that since the disorder is ('time') $z$-dependent, the system is of a different nature from the case where the interplay of thermalization and Anderson localization is studied [97, 98, 99, 100, 65, 102].

### 2.3.2 Vanishing correlations among modes

First of all, we study the correlations among the modes through the analysis of the evolution of the second-order moments in the 2×2 matrix $\langle a_p^* a_q^T \rangle (z)$. The impact of disorder is treated by making use of the Furutsu-Novikov theorem. This reveals that conservative disorder introduces an effective dissipation in the system. In this respect, we recall that, in principle, the laser beam excites strongly correlated modes at the fiber input. It is the effective dissipation due to disorder that breaks such a strong modal phase-correlation.

The theory developed in the Appendix A.1 reveals that in the regime discussed in Sec. 2.3.1, the correlations $\langle a_p^* a_q^T \rangle (z)$ among different modes ($p \neq q$), or within a single mode ($p = q$), both have the form

$$\langle a_p^* a_q^T \rangle (z) \sim \frac{\gamma (\mathbf{P}_p(\mathbf{a}(z)))^* a_q^T(z) - a_p^*(z) \mathbf{P}_q(\mathbf{a}(z))^T}{4\Delta\beta - i(\beta_p - \beta_q)},$$

Recalling that $L_d \ll L_{nl}$ we see that, at leading order, modal correlations vanish for propagation lengths larger than the nonlinear length, $\langle a_p^* a_q^T \rangle (z) \simeq 0$ for $z \gg L_{nl}$.

### 2.3.3 Closure of the moments equations

In the following, we derive the kinetic equation by following the wave turbulence perturbation expansion procedure, where linear dispersive effects dominate nonlinear effects $L_{lin} \ll L_{nl}$. Accordingly, an effective large separation of the linear and the nonlinear lengths scales takes place [40, 82, 42, 78, 41]. When combined with disorder effects, the modes exhibit random phases with quasi-Gaussian statistics, which allows one to achieve the closure of the infinite hierarchy of the moment equations. More precisely, because of the nonlinear character of the NLS equation, the evolution of the second-order moment
of the field depends on the fourth-order moment, while the equation for the fourth-order moment depends on the sixth-order moment, and so on. In this way, one obtains an infinite hierarchy of moment equations, in which the n-th order moment depends on the (n + 2)-th order moment of the field. This makes the equations impossible to solve, unless some method can be found to truncate the hierarchy. The closure of the infinite hierarchy of moment equations can be obtained in the weakly nonlinear regime by virtue of the Gaussian moment theorem. We remark in this respect that the key assumption underlying the wave turbulence approach is the existence of a random phase among the modes, rather than a genuine Gaussian statistics, as recently discussed in detail in Refs.[41, 130, 131]. We emphasize that in the present work the random phase of the modes is induced by the structural disorder of the medium, which dominates nonlinear effects ($L_d \ll L_{nl}$).

As discussed here above, the nondiagonal components of the 2×2 matrix $\langle a_p^* a_p^T \rangle (z)$ vanish. Accordingly, in the following we derive an equation governing the evolutions of the diagonal components $w_p(z) = \frac{1}{2} \langle |a_p(z)|^2 \rangle$. Starting from the modal NLS Eq.(3), we have

$$ \partial_z w_p = \frac{1}{3}\gamma \langle X_p^{(1)} \rangle + \frac{2}{3}\gamma \langle X_p^{(2)} \rangle, \quad (10) $$

$$ X_p^{(1)} = \operatorname{Im}\left\{ \sum_{l,m,n} S_{lmnp}^*(a_l^\dagger a_m^*) (a_n^T a_p) \right\}, \quad (11) $$

$$ X_p^{(2)} = \operatorname{Im}\left\{ \sum_{l,m,n} S_{lmnp}^*(a_n^T a_m^*) (a_l^\dagger a_p) \right\}. \quad (12) $$

The detailed derivation of the equations for the fourth order correlators $J_{lmnp}^{(1)}(z) = \langle (a_l^\dagger a_m^*)(a_n^T a_p) \rangle$ and $J_{lmnp}^{(2)}(z) = \langle (a_n^T a_m^*)(a_l^\dagger a_p) \rangle$ is given in the Appendix A.2. As already noticed, as a result of the Furutsu-Novikov theorem, conservative disorder introduces an effective dissipation in the system, so that the evolutions of the fourth-order moments have the form of a forced-damped oscillator equation:

$$ \partial_z J_{lmnp}^{(j)} = (-8\Delta\beta + i\Delta\omega_{lmnp}) J_{lmnp}^{(j)} + i\gamma \langle Y_{lmnp}^{(j)} \rangle \quad (13) $$

where $\Delta\omega_{lmnp} = \beta_l + \beta_m - \beta_n - \beta_p$ is the resonance frequency, while $\langle Y_{lmnp}^{(j)} \rangle$ denote the six-th order moments that have been split into products of second-order moments, by virtue of the factorizability property of statistical Gaussian fields (see the detailed expressions of $\langle Y_{lmnp}^{(j)} \rangle$ for $j = 1, 2$ in Eq.(A.8) and Eq.(A.15) in the Appendix A.2). The solution to Eq.(13) reads

$$ J_{lmnp}^{(j)}(z) = J_{lmnp}^{(j)}(0)G_{lmnp}(z) + \mathcal{I}_{lmnp}^{(j)}(z) \quad (14) $$

where the convolution integral reads

$$ \mathcal{I}_{lmnp}^{(j)}(z) = i\gamma \int_0^z \langle Y_{lmnp}^{(j)} \rangle (z-z') G_{lmnp}(z') dz' \quad (15) $$
with the Green function

$$G_{lmnp}(z) = H(z) \exp(i\Delta\omega_{lmnp}z - 8\Delta\beta z), \quad (16)$$

and $H(z)$ the Heaviside function. We now arrive at the key point of the analysis, in which one makes a distinction between the continuous and the discrete wave turbulence approaches.

### 2.3.4 Continuous wave turbulence

In general, an exchange of power among four different modes does not need to satisfy the exact resonance condition $\Delta\omega_{lmnp} = 0$. Indeed, it is sufficient for a mode quartet to verify the quasi-resonant condition

$$|\Delta\omega_{lmnp}| \lesssim 1/L_{kin}^{disor} \quad (17)$$

to provide a non-vanishing contribution to the convolution integral (15), where $L_{kin}^{disor}$ denotes the characteristic evolution length scale of the moments of the field $\langle Y_{lmnp}^{(j)} \rangle (z)$ (see Eq.(23) below). In the usual continuous wave turbulence approach, there is a large number of non-resonant mode quartets $\Delta\omega_{lmnp} \neq 0$ that verifies the quasi-resonant condition (17), and that contributes significantly to the evolution of the moments equations.

We anticipate that the continuous wave turbulence regime does not correspond to usual experiments of optical beam self-cleaning [85, 87], since we have typically $\beta_0 \sim 10^3 m^{-1}$, so that $\beta_0 \gg 1/L_{kin}^{disor}$. In other words, since any non-resonant mode quartet verifies $|\Delta\omega_{lmnp}| \ge \beta_0 \gg 1/L_{kin}^{disor}$, the system does not exhibit quasi-resonances, but solely exact resonances $\Delta\omega_{lmnp} = 0$. As a matter of fact, experiments of beam self-cleaning are described by the discrete wave turbulence approach that will be discussed in detail in the next Sec. 2.3.5. Here, for the sake of clarity, we first discuss the conventional continuous wave turbulence regime [41], which can be relevant for MMFs characterized by a huge number of modes ($M \gg 1$ such that $\beta_0 \ll \Delta\beta$). The discussion of the usual continuous regime is also important, in that it enlightens the effect of acceleration of thermalization, owing to the presence of structural disorder in the system.

Recalling that $L_d = 1/\Delta\beta \ll L_{nl}$, then the Green function decays on a length scale which is much smaller than the evolution length of $\langle Y_{lmnp}^{(j)} \rangle (z)$, so that the convolution integral can be approximated for $z \gg L_d$

$$J_{lmnp}^{(j)} \approx \gamma \langle Y_{lmnp}^{(j)} \rangle \frac{i8\Delta\beta - \Delta\omega_{lmnp}}{\Delta\omega_{lmnp}^2 + (8\Delta\beta)^2}. \quad (18)$$

By substitution of Eq.(18) in the fourth-order moments (14), one obtains the expressions of the averaged moments $\langle X_p^{(j)} \rangle$ given in (11-12). Collecting all terms in (10) gives the equation for the modal amplitudes $n_p(z) = 2w_p(z) = \langle |a_p(z)|^2 \rangle$. Since the MMF supports a large number of modes $M \gg 1$, we
consider the continuous limit where the discrete sums in (11-12) are replaced by continuous integrals. We obtain the continuous kinetic equation

$$ \partial_z \tilde{n}_{\mathbf{k}_4}(z) = \frac{4\gamma^2}{3\beta_0^6} \iiint d\mathbf{k}_{1,2,3} \frac{\overline{\Delta\beta}}{\Delta\tilde{\omega}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4}^2 + \overline{\Delta\beta}^2} \times |\tilde{S}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4}|^2 \tilde{M}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4} \\ + \frac{32\gamma^2}{9\beta_0^2} \int d\mathbf{k}_1 \frac{\overline{\Delta\beta}}{\Delta\tilde{\omega}_{\mathbf{k}_1\mathbf{k}_4}^2 + \overline{\Delta\beta}^2} |\tilde{s}_{\mathbf{k}_1\mathbf{k}_4}(\tilde{\mathbf{n}})|^2 (\tilde{n}_{\mathbf{k}_1} - \tilde{n}_{\mathbf{k}_4}) \quad (19) $$

where $d\mathbf{k}_{1,2,3} = d\mathbf{k}_1 d\mathbf{k}_2 d\mathbf{k}_3$, $\overline{\Delta\beta} = 8\Delta\beta$ and

$$ \tilde{s}_{\mathbf{k}_1 \mathbf{k}_4}(\tilde{\mathbf{n}}) = \frac{1}{\beta_0^2} \int d\mathbf{k}' \tilde{s}_{\mathbf{k}_1 \mathbf{k}' \mathbf{k}' \mathbf{k}_4} \tilde{n}_{\mathbf{k}'} \quad (20) $$

The functions with a tilde refer to the natural continuum extension of the corresponding discrete functions, i.e., $\tilde{n}_{\mathbf{k}}(z) = n_{[\mathbf{k}/\beta_0]}(z)$, $\tilde{\beta}_{\mathbf{k}} = \beta_{[\mathbf{k}/\beta_0]}$, $\tilde{S}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4} = S_{[\mathbf{k}_1/\beta_0][\mathbf{k}_2/\beta_0][\mathbf{k}_3/\beta_0][\mathbf{k}_4/\beta_0]}$ and so on, where $[x]$ denotes the integer part of $x$. With these notations we have $\tilde{M}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4} = \tilde{n}_{\mathbf{k}_1}\tilde{n}_{\mathbf{k}_2}\tilde{n}_{\mathbf{k}_3}\tilde{n}_{\mathbf{k}_4}(\tilde{n}_{\mathbf{k}_4}^{-1} + \tilde{n}_{\mathbf{k}_3}^{-1} - \tilde{n}_{\mathbf{k}_1}^{-1} - \tilde{n}_{\mathbf{k}_2}^{-1})$, $\Delta\tilde{\omega}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4} = \tilde{\beta}_{\mathbf{k}_1} + \tilde{\beta}_{\mathbf{k}_2} - \tilde{\beta}_{\mathbf{k}_3} - \tilde{\beta}_{\mathbf{k}_4}$, $\Delta\tilde{\omega}_{\mathbf{k}_1\mathbf{k}_2} = \tilde{\beta}_{\mathbf{k}_1} - \tilde{\beta}_{\mathbf{k}_2}$, and $\tilde{\beta}_{\mathbf{k}} = \kappa_x + \kappa_y + \beta_0$, with $\kappa = \beta_0(p_x, p_y)$ [11].

The main novelty of the kinetic Eq.(20) with respect to the case without any structural disorder [11], is that the mechanism of disorder-induced dissipation introduces a finite bandwidth into the four-wave resonances among the modes. Accordingly, instead of the Dirac $\delta$-distribution that guarantees energy conservation at each four-wave interaction, here the kinetic equation involves a Lorentzian distribution. We remark that this aspect was already discussed in different circumstances [77], in particular in recent work [71] dealing with random fiber lasers in the presence of gain and losses. Here, the originality is that the finite bandwidth of the Lorentzian distribution and the associated effective dissipation $\Delta\beta$ originate from the conservative structural disorder of the material.

Taking the formal limit $\Delta\beta \to 0$: $\overline{\Delta\beta}/(\Delta\tilde{\omega}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4}^2 + \overline{\Delta\beta}^2) \to \pi\delta(\Delta\tilde{\omega}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4})$, we recover a continuous kinetic equation with a form similar to that derived in [11] in the absence of disorder ($\Delta\beta = 0$) and in the absence of polarization effects (scalar limit $a_p \to a_{p,x}$). However, the above limit $\Delta\beta \to 0$ is not physically relevant here, since we have assumed $L_d = 1/\Delta\beta \ll L_{nl}$ in order to derive Eq.(13). Here, we are considering the regime $\beta_0 \ll \Delta\beta$, so that it is the opposite limit that is physically relevant $\overline{\Delta\beta}/(\Delta\tilde{\omega}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4}^2 + \overline{\Delta\beta}^2) \to 1/\overline{\Delta\beta}$, and the kinetic equation takes the reduced form:

$$ \partial_z \tilde{n}_{\mathbf{k}_4}(z) = \frac{\gamma^2}{6\Delta\beta\beta_0^6} \iiint d\mathbf{k}_{1,2,3} |\tilde{S}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4}|^2 \tilde{M}_{\mathbf{k}_1\mathbf{k}_2\mathbf{k}_3\mathbf{k}_4} \\ + \frac{4\gamma^2}{9\Delta\beta\beta_0^2} \int d\mathbf{k}_1 |\tilde{s}_{\mathbf{k}_1\mathbf{k}_4}(\tilde{\mathbf{n}})|^2 (\tilde{n}_{\mathbf{k}_1} - \tilde{n}_{\mathbf{k}_4}) \quad (21) $$

The second term in the right-hand side of (21) enforces an isotropic distribution of mode occupancies $\tilde{n}_{\mathbf{k}}(z)$ for degenerate modes, while the first term enforces
the mode occupancies to reach the most disordered equilibrium distribution. The kinetic equation (21) conserves the power, $N = \beta_0^{-2} \int \tilde{n}_{\mathbf{k}} d\mathbf{k}$, and exhibits a H-theorem of entropy growth, $\partial_z S^{\text{neq}}(z) \ge 0$, where the nonequilibrium entropy is defined by $S^{\text{neq}}(z) = \beta_0^{-2} \int \log(\tilde{n}_{\mathbf{k}}) d\mathbf{k}$. However, at variance with the conventional wave turbulence kinetic equation, the kinetic Eq.(21) does not conserve the energy, $E = \beta_0^{-2} \int \tilde{\beta}_{\mathbf{k}} \tilde{n}_{\mathbf{k}} d\mathbf{k}$. Then the equilibrium distribution maximizing the entropy, given the constraint of the conservation of power $N$, is given by the uniform distribution

$$ \tilde{n}_{\mathbf{k}}^{\text{eq}} = \text{const.} \qquad (22) $$

This equilibrium state denotes an equipartition of ('particles') power among all modes. Recalling that $\beta_0 \ll \Delta\beta$, the Lorentzian distribution in the kinetic Eq.(19) is dominated by disorder, so that the length scale characterizing the rate of thermalization toward the equilibrium distribution (22) is given by

$$ L_{kin}^{\text{disor}} \sim \Delta\beta L_{nl}^2 / \bar{S}_{lmnp}^2, \qquad (23) $$

where $\bar{S}_{lmnp}^2$ denotes the average square of the tensor $S_{lmnp}$ involving non-trivial resonances among nondegenerate modes.

The main conclusion is that the kinetic Eq.(20) does not describe a process of condensation characterized by a macroscopic occupation of the fundamental mode. This indicates that weak disorder should prevent an effect of beam self-cleaning in MMFs possessing a huge number of modes.

### 2.3.5 Discrete wave turbulence

We have seen that, in the continuous wave turbulence regime, quasi-resonances verifying $\Delta\omega_{lmnp} \lesssim 1/L_{kin}^{\text{disor}}$ contribute to the convolution integral (15). At variance with the continuous regime, in the discrete case the non-vanishing minimum value of $\Delta\omega_{lmnp}$ is such that

$$ \min(|\Delta\omega_{lmnp}|) = \beta_0 \gg 1/L_{kin}^{\text{disor}}. \qquad (24) $$

Accordingly, only exact resonances $\Delta\omega_{lmnp} = 0$ contribute to the integral (15), while non-resonant mode quartets lead to a vanishing integral [41]. Note that this procedure can also be justified by a homogenization approach, as reported in [132] in the presence of a non-conservative disorder accounting for gain and losses in the system. As discussed above, usual experiments of optical beam self-cleaning are described by the discrete wave turbulence approach, since $\beta_0^{-1} \sim 10^{-3} \text{m}$ and the above condition is well verified in the experiments [85, 87].

In the discrete wave turbulence regime we need to consider separately the cases of resonant and non-resonant mode interactions. For mode quartets verifying $\Delta\omega_{lmnp} = 0$, the Green function (16) decays on a length scale which is much smaller than the evolution of $n_p(z)$, because $L_d = (\Delta\beta)^{-1} \ll L_{kin}$, so that (15) can be approximated by $J_{lmnp}^{(j)}(z) \approx \frac{i\gamma}{8\Delta\beta} \langle Y_{lmnp}^{(j)} \rangle(z)$. For $\Delta\omega_{lmnp} \ne 0$, the
Green function oscillates rapidly over a length scale smaller than $L_{lin} = \beta_0^{-1} \ll L_{kin}$. Such a rapid phase rotation, combined with the fast decay of the Green function over a length $L_d \ll L_{nl}$, leads to a vanishing convolution integral in (15). As a result, the fourth-order moment can be written in the form

$$J_{lmnp}^{(j)}(z) \approx \frac{i\gamma}{8\Delta\beta} \langle Y_{lmnp}^{(j)} \rangle (z) \delta^K(\Delta\omega_{lmnp}), \quad (25)$$

where $\delta^K(\Delta\omega_{lmnp}) = 1$ if $\Delta\omega_{lmnp} = 0$, and zero otherwise. Note that the discrete regime discussed here does not exactly correspond to the discrete regimes associated with finite size effects in homogeneous wave turbulence [41] – here the system is non-homogeneous ($V(\mathbf{r}) \neq \text{const}$) and the resonance for the momentum reflected by the tensor $S_{lmnp}$ is not as rigid as the usual one involving the Dirac $\delta$-distribution in homogeneous turbulence.

We provide a detailed computation of the fourth-order moments $J_{lmnp}^{(j)}(z)$ and corresponding moments $\langle X_p^{(j)} \rangle$ defined in (11-12) in the Appendix A.2, see Eq.(A.12) and (A.17). Then collecting all terms in (10) gives the discrete kinetic equation for the modal amplitudes $n_p(z) = 2w_p(z)$

$$\partial_z n_p(z) = \frac{\gamma^2}{6\Delta\beta} \sum_{l,m,n} |S_{lmnp}|^2 \delta^K(\Delta\omega_{lmnp}) M_{lmnp}(\mathbf{n}) \\ + \frac{4\gamma^2}{9\Delta\beta} \sum_l |s_{lp}(\mathbf{n})|^2 \delta^K(\Delta\omega_{lp})(n_l - n_p), \quad (26)$$

with $s_{lp}(\mathbf{n}) = \sum_{m'} S_{lm'm'p} n_{m'}$, and $M_{lmnp}(\mathbf{n}) = n_l n_m n_p + n_l n_m n_n - n_n n_p n_m - n_n n_p n_l$, with 'n$_m$' for 'n$_m$'($z$), $\Delta\omega_{lp} = \beta_l - \beta_p$. According to the kinetic equation (26), the length scale characterizing the rate of thermalization is the same as that obtained in the continuous wave turbulence regime (23), namely $L_{kin}^{disor} \sim \Delta\beta L_{nl}^2 / \bar{S}_{lmnp}^2$.

Aside from its discrete form, the kinetic Eq.(26) has a structure which is analogous to the conventional wave turbulence kinetic equation, so that it describes a process of wave condensation that occurs irrespective of the sign of the nonlinearity $\gamma$ (see the factor $\gamma^2$ in Eq.(26))[57, 41]. We refer the reader to Refs.[57, 133, 41] for details on condensation in the homogeneous case ($V(\mathbf{r}) = 0$), and to [11, 134] for the non-homogeneous case in a waveguide potential ($V(\mathbf{r}) \neq 0$). It is important to note that, in contrast to the kinetic Eq.(21) derived in the continuous wave turbulence regime, here Eq.(26) also conserves the energy $E = \sum_p \beta_p n_p(z)$ – despite the presence of the dissipation effect $\Delta\beta$. The reason for this is that only exact resonances contribute to the discrete turbulence regime, so that the discrete kinetic equation is not affected by the dissipation-induced resonance broadening which inhibits energy conservation in the continuous turbulence regime. The kinetic Eq.(26) also conserves the 'number of particles' $N = \sum_p n_p(z)$ and exhibits a H-theorem of entropy growth

$$\partial_z S^{neq} \geq 0, \quad (27)$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Scaling of acceleration of thermalization with mode-decorrelated disorder: Numerical simulations of the NLS Eq.(3) showing the evolutions of the fundamental mode $n_0(z)$, for different amounts of disorder $\Delta\beta$. The dashed lines show the corresponding simulations of the discrete kinetic Eq.(26), starting from the same initial condition as the NLS simulations. Parameters are: $\Delta\beta \simeq 2.6 \text{ m}^{-1}$ ($2\pi/\sigma_\beta = 2.1 \text{ m}$, $l_\beta = 30 \text{ cm}$) blue (dash-dotted); $\Delta\beta \simeq 10.5 \text{ m}^{-1}$, ($2\pi/\sigma_\beta = 26 \text{ cm}$, $l_\beta = 1.88 \text{ cm}$) green (dotted); $\Delta\beta \simeq 42 \text{ m}^{-1}$, ($2\pi/\sigma_\beta = 6.6 \text{ cm}$, $l_\beta = 0.47 \text{ cm}$) red (solid); while the power is $N = 47.5 \text{ kW}$ ($M = 120$ modes, fiber core radius $R = 26 \text{ µm}$). The curves eventually relax to the common theoretical equilibrium value $n_0^{eq}/N \simeq 0.2$ (dashed black line) with different rates, confirming the scaling of acceleration of condensation predicted by the theory in Eq.(30), without using adjustable parameters. Source: From Ref.[25].</figcaption></figure>

for the nonequilibrium entropy

$$S^{neq}(z) = \sum_{p} \log(n_{p}(z)). \qquad (28)$$

The growth of entropy saturates once equilibrium is reached, as discussed in a phenomenological way in section 2.1. Accordingly, the kinetic Eq.(26) describes an irreversible evolution to the (maximum entropy) thermodynamic Rayleigh-Jeans equilibrium

$$n_{p}^{eq} = T/(\beta_{p} - \mu). \qquad (29)$$

Note that this equilibrium distribution coincides with the equilibrium (2) discussed in sec.2.1. Indeed, the equilibrium distribution (29) does not depend on disorder. However, disorder impacts the rate of thermalization to this equilibrium, a property that is discussed in the next section.

### 2.3.6 Acceleration of thermalization mediated by disorder

The length scale characterizing the rate of thermalization in the absence of structural disorder is obtained from the continuous kinetic equation that was derived in Ref.[11], $L_{kin}^{ord} \sim \beta_0 L_{nl}^2 / \bar{S}_{lmnp}^2$. A similar scaling behavior is expected in the discrete wave turbulence regime, see [25]. Hence, in regimes of beam cleaning experiments where $\beta_0 \gg \Delta\beta$, the weak disorder leads to a significant acceleration of the rate of thermalization and condensation

$$L_{kin}^{disor} / L_{kin}^{ord} \sim \Delta\beta/\beta_0. \qquad (30)$$
Considering typical values of $\beta_0 \sim 10^3 m^{-1}$ and $L_d = 1/\Delta\beta$ larger than tens of centimeters, we see that we always have $L_{kin}^{disor}/L_{kin}^{ord} \ll 1$. Accordingly, weak disorder is responsible for a significant acceleration of the rate of thermalization and condensation [103]. Note that the presence of a perturbation on the dispersion relation can modify the above scaling, a feature that will be discussed later on.

We have confirmed the scaling (30) by performing numerical simulations of the NLS Eq.(3) for different amounts of disorder $\Delta\beta$. The results are reported in Fig. 3, and confirm the scaling of the rate of acceleration of thermalization and condensation as they are predicted by the theory. More precisely, a remarkable quantitative agreement has been obtained between the simulations of the NLS Eq.(3) and those of the discrete kinetic Eq.(26), without using any adjustable parameter. Note that the equilibrium value of the condensate fraction $n_0^{eq}/N$ vs. the energy $E$ is obtained from the analysis of the RJ equilibrium distribution [57, 11, 41]. The condensation curve $n_0^{eq}/N$ vs. the energy $E$ is reported in explicit form in the Appendix B, and its experimental study is discussed in section 3.3.1. Also note that we have deliberately chosen a small value of the condensate fraction $n_0/N \simeq 0.2$, so as to avoid large deviations from Gaussianity for the fundamental mode – though the theory has been validated even for large condensate fractions in [103] (also see Fig. 4 here below). We recall here the recent works showing that a key assumption of the wave turbulence approach is the existence of a random phase among the modes, rather than a genuine Gaussian statistics [130, 131].

### 2.3.7 Light thermalization resulting from turbulence cascades

Spatial beam self-cleaning in MMFs can also be understood in the context of turbulence cascades associated to conserved quantities, as initially explored in the broader framework of hydrodynamic turbulence in Ref.[136]. To properly contextualize this problem within the wave turbulence theory [40, 41, 42], let us consider a generic dissipative system that is driven far from equilibrium by an external source. A typical illustrative example of forced physical system can be the excitation of hydrodynamic surface waves by the wind in oceans. In general, the frequency-scales of forcing and damping differ significantly. The nonlinear interaction leads to a spectral redistribution, which eventually relaxes to a nonequilibrium stationary state characterized by a non-vanishing flux of the conserved quantities, typically the power (wave action), and the kinetic energy. The nonequilibrium stationary solutions of the wave turbulence kinetic equation were originally obtained by Zakharov in the mid-1960s [137, 138], and they can be considered as the analogue of the spectra of hydrodynamic turbulence proposed in 1941 by Kolmogorov [139].

Although nonequilibrium stationary spectra of turbulence are sustained by a continuous forcing and damping, it is important to remark that they can also emerge in the transient evolution of a purely conservative system, prior to reaching the RJ equilibrium distribution, see e.g. Refs.[59, 60, 140]. It is from this perspective that the irreversible process of light thermalization in MMFs has been addressed in Ref.[135]. In this article, numerical simulations of the
<figure><img src="image_5.png" /><figcaption>Figure 4: Energy and power flows underlying light thermalization: Numerical simulation of the modal NLS Eq.(3) (a)-(b), and wave turbulence kinetic Eq.(26) (c)-(d): Evolutions of the modal distribution of the power $\tilde{n}_g$ (a)-(c), and energy $\tilde{\mathcal{E}}_g$ (b)-(d), for $g_{max}$ =15 groups of non-degenerate modes of the parabolic MMF (an average has been taken among the degenerate modes). The dashed black lines in (a) and (c) denote the RJ power-law $\tilde{n}_g^{eq} \sim 1/g$. The thermalization to the RJ equilibrium is characterized by an energy flow toward the higher-order modes, and a power flow toward the fundamental mode and the higher-order modes. Evolutions of $n_0(z)$ (e) and $n_4(z)$ (f) obtained from simulation of the NLS Eq.(3) (red line) and the kinetic Eq.(26) (dashed blue): The modal components thermalize to the theoretical equilibrium value predicted by the RJ theory (the dashed black line denotes $n_0^{eq}/N = 0.6$). Parameters: $N = 47.5\,\text{kW}$, $l_\beta = 0.019\,\text{m}$, $2\pi/\sigma_\beta = 0.26\,\text{m}$, there is no average over the realizations for the NLS simulation. Source: From Ref.[135].</figcaption></figure>
NLS Eq.(3), and the corresponding kinetic Eq.(26), have been realized by start-
ing from the same initial mode distribution (except that independent random
phases among the modes were considered for NLS simulations). The results
are reported in Fig. 4. During light propagation in the MMF, the spectrum
for the power distribution $n_p$ essentially flows toward the fundamental mode
(inverse cascade), while a small fraction of $n_p$ flows toward the higher-order
modes. For convenience, we have reported in Fig. 4 the average power $\tilde{n}_g$ within
each group of degenerate modes, where $g = 0, \dots, g_{max} - 1$ indexes the mode
group (in Fig. 4 $g_{max} = 15$ for a total $M = g_{max}(g_{max} + 1)/2 = 120$ modes).
The RJ power-law $\tilde{n}_g \sim 1/g$ is verified by the simulation of the kinetic and
NLS equations. Note that, due to the extensive computation times required,
averaging over the realizations of NLS simulations was not feasible, resulting
in the spectrum's noisy structure of $\tilde{n}_g$ in Fig. 4(a)-(b).

The results for the power mode distribution are corroborated by those of the
modal distribution of the energy, which exhibits a flow toward the higher-order
modes (direct cascade). In this example, we considered a relatively small value
of the conserved energy *E*, which is below the critical value of the transition to
condensation *E*<sub>crit</sub>. In the condensed state, the chemical potential tends to the
lowest energy level, μ → β<sub>0</sub><sup>-</sup> (see Appendix B), so that the waves that started
from an initial state with an excess energy in the low-energy modes, eventually
tend to an equilibrium state displaying an energy equipartition among the
modes ε<sub>p</sub> = (β<sub>p</sub> − β<sub>0</sub>)n<sub>p</sub> ∼ *T* [or equivalently Ė<sub>g</sub> = β<sub>0</sub>g*tilde{n}*<sub>*g*</sub> ∼ *T*], as illustrated
in Fig. 4(b)-(d). Then RJ thermalization is characterized by a macroscopic
population of the fundamental mode, as illustrated in Fig. 4(e), where the
condensate fraction relaxes toward the theoretical equilibrium value n<sub>0</sub><sup>eq</sup>/N ≃
0.6. Note that the good agreement between the simulations of the NLS and
kinetic equations in Fig. 4 is obtained without using adjustable parameters.

One may question whether the above power and energy flows can be de-
scribed theoretically by means of nonequilibrium stationary solutions of the
kinetic equation – the so-called Zakharov-Kolmogorov spectra of turbulence
[40]. In this regard, it is crucial to note that the kinetic Eq.(26) considered
here, differs from the conventional wave turbulence kinetic equation in two
key aspects: (i) It involves the tensor |S<sub>plmn</sub>|<sup>2</sup> instead of the Dirac δ–function
over the wave-vectors, because the transverse trapping potential inherent to
the MMF waveguide breaks the conservation of the transverse momentum. (ii)
The kinetic Eq.(26) is discrete in frequencies. This latter property does not
allow the application of the standard procedure based on the Zakharov con-
formal transformation to derive nonequilibrium stationary solutions featured
by a non-vanishing flux of the conserved power and energy. These remarks
appear consistent with the results of the numerical simulations of the NLS and
kinetic equations, which do not evidence the formation of a nonequilibrium
power-law spectrum in the transient evolution that precedes the formation of
the equilibrium RJ spectrum, ñ<sub>g</sub><sup>eq</sup> ∼ 1/g.

We finally note that an experimental signature of the direct and inverse
turbulence cascades underlying RJ thermalization in MMFs has been recently
<figure><img src="image_6.png" /><figcaption>Figure 5: Disorder-induced beam cleaning: Numerical simulations of the NLS Eq.(3) showing the evolution of the relative variance of intensity fluctuations $\sigma_I(z)$ from Eq.(31) (red circles), and of the condensate fraction $n_0(z)/N$ (blue line), in the presence of disorder (a), and in the absence of disorder (b). The intensity patterns $I(r)$ inside the fiber core (circles) are shown at $z = 0$ and $z = 15$ m. The horizontal dashed black line in (a) shows the condensate fraction at thermal equilibrium, $n_0^{eq}/N \simeq 0.68$. The disorder induces a stable beam cleaning condensation, characterized by a significant reduction of intensity fluctuations down to $\sigma_I \simeq 0.1$ (a), which is in contrast with the case without disorder, where the variance of intensity fluctuations remain almost constant throughout propagation with $\sigma_I \simeq 0.9$. Parameters are: $\Delta\beta = 2.6$ m ($2\pi/\sigma_\beta = 2.13$ m, $l_\beta = 30$ cm), $M = 120$ modes, fiber core radius $R = 26$ µm, $N = 47.5$ kW. Source: From Ref.[25].</figcaption></figure>

reported in Ref.[141]. More precisely, it was shown that light thermalization is characterized by a flow of the energy toward the higher-order modes, and a bi-directional redistribution of power from intermediate modes toward the fundamental mode and the higher-order modes.

### 2.3.8 Disorder-induced beam cleaning

The amount of 'beam cleaning' can be quantified through the analysis of the fluctuations of the intensity during the propagation in the MMF. To this aim, we consider the relative variance of intensity fluctuations relevant for a spatially non-homogeneous incoherent beam

$$ \sigma_I^2(z) = \frac{\int \langle I^2(\mathbf{r}, z) \rangle - \langle I(\mathbf{r}, z) \rangle^2 d\mathbf{r}}{\int \langle I(\mathbf{r}, z) \rangle^2 d\mathbf{r}}, \quad (31) $$
where the intensity is $I(\mathbf{r}, z) = |\psi|^2(\mathbf{r}, z)$. Note that for a beam with Gaussian statistics we have $\sigma_I^2 = 1$. We report in Fig. 5 the evolutions of the variance of intensity fluctuations in the presence and the absence of disorder, respectively. The brackets $\langle \cdot \rangle$ in Eq.(31) refers to an averaging over the propagation length $\Delta z = 10$ mm, which is larger than the mode-beating length scale $\sim \beta_0^{-1}(\approx 0.2$mm in the example of Fig. 5). The simulations in Fig. 5 clearly show that the presence of disorder induces a rapid condensation process, which in turn leads to a significant reduction of the relative standard deviation of intensity fluctuations $\sigma_I \approx 0.1$. Conversely, in the absence of disorder the multimode beam can exhibit an enhanced brightness at some propagation lengths (see the intensity pattern at $z = 15$ m in Fig. 5(b)), however its oscillatory multimode nature prevents a stable beam-cleaning propagation, as evidenced by the relative standard deviation of intensity fluctuations that only slowly decreases below $\sigma_I \approx 0.9$.

### 2.3.9 Correlated and partially correlated disorder

In the previous section 2.3.1 we have considered a model of disorder that can be termed 'mode-decorrelated', in the sense that each individual mode of the MMF experiences a different noise, i.e., the functions $\nu_{p,j}(z)$ in (6) are independent of each other. Although this approach can be considered to be justified in different circumstances [116, 117, 142], one may question its validity for the case of an MMF exhibiting a large number of modes. In the following we address this question by considering two different models of disorder, namely the fully mode-correlated disorder model, and the partially mode-correlated disorder model.

**Mode-correlated noise:** We first consider the fully mode-correlated model that can be considered as the opposite limit of the 'mode-decorrelated' one, in the sense that all modes experience the same noise (more precisely, the same realization of the noise). In this limit, the $2 \times 2$ matrices describing the modal noise in (6) reduce to $D_p = D$ with $D = \sum_{j=0}^{3} \nu_j \sigma_j$.

The theory developed above for the decorrelated model of disorder can be extended to this fully correlated model. The theory reveals a remarkable result in this case, namely, that the equations for the fourth-order moments $J_{lmnp}^{(j)}(z)$ do not exhibit an effective damping, i.e., $\Delta\beta = 0$, see Appendix A.3. This is in marked contrast with the previously discussed decorrelated model of disorder, see Eqs.(13). This result has a major consequence: the fully mode-correlated noise does not modify the regularization of wave resonances, and the system recovers a dynamics which is analogous to that obtained in the absence of any disorder. We illustrate this in Fig. 6, that reports the evolution of the modal components $n_p(z)$ obtained by simulation of the modal NLS Eq.(3) in the presence of a mode-correlated disorder. The initial condition is coherent (all modes are correlated with each other), and we see in Fig. 6 that the low-order modes recover an oscillatory dynamics reflecting the presence of strong phase-correlations, as for the coherent regime discussed in the absence of disorder, see Fig. 2. This coherent dynamics is consistent with the intuitive idea that, in
<figure><img src="image_7.png" /><figcaption>Figure 6: Mode-correlated disorder: Numerical simulations of the NLS Eq.(3) showing the evolutions of the modal components $n_p(z)$, starting from a coherent initial condition: fundamental mode $p = 0$ (red dashed), $p = 1$ (dark blue solid), $p = 2$ (blue solid), $p = 3$ (light blue solid), $p = 4$ (cyan solid), $p = 5$ (light green solid), $p = 6$ (green solid), $p = 7$ (yellow solid), $p = 8$ (orange solid). The mode-correlated disorder does not introduce an effective dissipation in the system ($\Delta\beta = 0$): phase-correlations among the low-order modes are preserved, and lead to an oscillatory dynamics similar to that in the absence of disorder (see Fig. 2). At complete thermal equilibrium the system would reach the condensate fraction $n_0^{eq}/N \simeq 0.68$. Parameters are: $\Delta\beta \simeq 2.6 \text{ m}^{-1}$ ($2\pi/\sigma_\beta = 2.1 \text{ m}$, $l_\beta = 30 \text{ cm}$), the power is $N = 47.5 \text{ kW}$, $M = 120$ modes, fiber core radius $R = 26 \text{ µm}$. Source: From Ref.[25].</figcaption></figure>

the absence of an effective dissipation ($\Delta\beta = 0$), the phase-correlations are not forgotten during the evolution.

*Partially mode-correlated noise:* We have seen that the fully mode-correlated disorder does not introduce an effective dissipation ($\Delta\beta = 0$), and thus leads to a coherent dynamics which is analogous to that obtained in the absence of disorder. We have thus considered a 'partially correlated' model of disorder, in which modes that belong to different groups of degenerate modes experience a decorrelated noise, while degenerate modes of the same group experience the same noise. This artificial model for MMFs can be considered as an intermediate model between the two limits of correlated and decorrelated models. It is mathematically tractable, and it illustrates the conjecture that the discrete kinetic Eq.(26) is robust as soon as disorder is not fully mode-correlated.

The theory developed above for the model of decorrelated disorder has been extended to address the partially correlated model, see Appendix A.4. The theory reveals that (second-order) correlations among non-degenerate modes are vanishing small and can be neglected, as it was shown for the model of decorrelated disorder. However, the computation of the fourth-order moments $J_{lmnp}^{(j)}$ is more delicate, because the mode-correlated noise introduces more terms in the calculation of the equations for the moments. Different results for the fourth-order moments are obtained, which depend on the specific modes involved in the moments. Almost all of the fourth-order moments satisfy an evolution equation with a dissipation that is proportional to $\Delta\beta$. This result is analogous to that obtained for the model of decorrelated disorder considered above, though the coefficients in front of $\Delta\beta$ are different, and their values depend on the specific modes involved in the fourth-order moment. In
<figure><img src="image_8.png" /><figcaption>Figure 7: Scaling of acceleration of thermalization with a partially mode-correlated disorder: Numerical simulations of the NLS Eq.(3) showing the evolutions of the fundamental mode $n_0(z)$, for different amounts of disorder $\Delta\beta$. The dashed lines show the corresponding simulations of the discrete kinetic Eq.(26), starting from the same initial condition as the NLS simulations. Parameters are: $\Delta\beta \simeq 2.6 \text{ m}^{-1}$ ($2\pi/\sigma_\beta = 2.1$ m, $l_\beta = 30$ cm) blue (dashdotted); $\Delta\beta \simeq 10.5 \text{ m}^{-1}$, ($2\pi/\sigma_\beta = 26$ cm, $l_\beta = 1.88$ cm) green (dotted); $\Delta\beta \simeq 42 \text{ m}^{-1}$, ($2\pi/\sigma_\beta = 6.6$ cm, $l_\beta = 0.47$ cm) red (solid); while the power is $N = 47.5$ kW ($M = 120$ modes, fiber radius $R = 26$ $\mu$m). The curves eventually relax to the common theoretical equilibrium value $n_0^{eq}/N$ with different rates, and confirm the theoretical scaling of acceleration of condensation in Eq.(30) for a partially mode-correlated disorder. Source: From Ref.[25].</figcaption></figure>

addition, there are also particular cases where the fourth-order moments do not exhibit any dissipation. Such special cases do not contribute to the fast thermalization process described by the effective dissipation $\Delta\beta$, but instead they induce a reversible exchange of power within a group of degenerate modes, see Appendix A.4.

To summarize, the theoretical developments reported in Appendix A.4 allow to infer that the kinetic equation is still of the form given by the discrete kinetic Eq.(26), and that the scaling of the rate of thermalization is still given by $L_{kin}^{disor} \sim \Delta\beta L_{nl}^2 / \bar{S}_{lmnp}^2$ in Eq.(23). This scaling of the rate of thermalization has been confirmed by the numerical simulations of the NLS Eq.(3). The results are reported in Fig. 7 for the same parameters of disorder as those considered in Fig. 3, except that a model of partial mode-correlated disorder has been considered. The agreement with the discrete kinetic Eq.(26) confirms the scaling (30) of the rate of acceleration of thermalization for the 'partially correlated' model of disorder.

### 2.3.10 Corrections of the dispersion relation

The wave turbulence approach developed in this section is based on the NLS equation with an ideal parabolic potential and the (linear) dispersion relation $\beta_p = \beta_0(p_x + p_y + 1)$. Several factors introduce perturbations to this expression of the dispersion relation, which can be written in the form $\tilde{\beta}_p = \beta_p + b_p$, where the perturbation $b_p$ is a function of $(p_x, p_y)$ with $b_0/\beta_0 \ll 1$. An example of perturbation is the leading order correction due to angular dispersion effects in the Helmholtz equation, $b_p = (\beta_0^2/(2k_0n_0))(1 + p_x + p_y)^2$, which typically can give
$b_0/\beta_0 \sim 5 \times 10^{-4}$ with usual beam-cleaning experimental parameters. As a matter of fact, such a correction to the dispersion relation was considered in [136] to study cascaded energy transfers among the fiber modes. Another example is provided by the well-known fact that a GRIN MMF usually exhibits deviations from the ideal parabolic shape, i.e., $V(\mathbf{r}) \sim |\mathbf{r}|^\nu$ with an exponent that deviates from $\nu = 2$. The general expression of the eigenvalue is rather complicated: it is of the form $\tilde{\beta}_p \propto (1+p_x+p_y)^{2\nu/(\nu+2)}$ [143]. Considering a deviation of a few percents from $\nu = 2$ [144], one has $b_0/\beta_0 \sim 5 \times 10^{-3}$ with typical parameters of beam cleaning experiments. In addition, we can notice that the truncation of the parabolic refractive index profile due to the fiber cladding can introduce perturbations of the higher-order propagation constants. It is also important to note that the standard deviation of the fluctuations of the structural disorder of the MMF due to imperfections and external perturbations (term $D_p a_p$ in the modal NLS Eq.(3)) may be of the same order as the correction of the dispersion relation.

Let us discuss the possible impact of perturbations of the dispersion relations. In this respect, resonances that are exact at leading order ($\Delta\omega_{lmnp} = \beta_l + \beta_m - \beta_n - \beta_p = 0$) exhibit a residual non-resonant contribution, i.e., $\Delta\tilde{\omega}_{lmnp} = \tilde{\beta}_l + \tilde{\beta}_m - \tilde{\beta}_n - \tilde{\beta}_p = \Delta b_{lmnp}$ with $\Delta b_{lmnp} = b_l + b_m - b_n - b_p$. The kinetic equation accounting for the correction on the dispersion relation can be derived with the above assumption $L_d = 1/\Delta\beta \ll L_{nl} < L_{disor}^{kin}$. The result of the convolution integral Eq.(15) is approximated by

$$J_{lmnp}^{(j)} \approx \gamma \langle Y_{lmnp}^{(j)} \rangle \frac{i8\Delta\beta - \Delta b_{lmnp}}{\Delta b_{lmnp}^2 + (8\Delta\beta)^2} \delta^K(\Delta\omega_{lmnp}). \quad (32)$$

Proceeding as in section 2.3.5, we obtain the discrete kinetic equation

$$\begin{aligned} \partial_z n_p(z) = & \frac{4\gamma^2 \overline{\Delta\beta}}{3} \sum_{l,m,n} \frac{\delta^K(\Delta\omega_{lmnp})}{\Delta b_{lmnp}^2 + \overline{\Delta\beta}^2} |S_{lmnp}|^2 M_{lmnp}(\mathbf{n}) \\ & + \frac{32\gamma^2 \overline{\Delta\beta}}{9} \sum_l \frac{\delta^K(\Delta\omega_{lp})}{\Delta b_{lp}^2 + \overline{\Delta\beta}^2} |s_{lp}(\mathbf{n})|^2 (n_l - n_p) \end{aligned} \quad (33)$$

where we recall that $s_{lp}(\mathbf{n}) = \sum_{m'} S_{lm'm'n'm'} n_{m'}$, $M_{lmnp}(\mathbf{n}) = n_l n_m n_p + n_l n_m n_n - n_n n_p n_m - n_n n_p n_l$, with $\Delta b_{lp} = b_l - b_p$ and $\overline{\Delta\beta} = 8\Delta\beta$. As already commented through Eq.(19), the Lorentzian distribution reflects the finite bandwidth of the four-wave resonances owing to the effective dissipation $\Delta\beta$. Accordingly, the kinetic Eq.(33) conserves $E = \sum_p \beta_p n_p(z)$, but not $\tilde{E} = \sum_p \tilde{\beta}_p n_p(z)$. In addition, the conservation of power $N = \sum_p n_p(z)$ and the H-theorem of entropy growth for $S(z) = \sum_p \log(n_p(z))$ describe a relaxation to $n_p^{eq} = T/(\beta_p - \mu)$, i.e., the same equilibrium as in the absence of the correction on the dispersion relation ($b_p = 0$). In summary, in the regime $\Delta\beta \gg b_0$, the Lorentzian distribution can be simplified $\overline{\Delta\beta}/(\Delta b_{lmnp}^2 + \overline{\Delta\beta}^2) \to 1/\overline{\Delta\beta}$, and the kinetic Eq.(33) exactly recovers the previous kinetic Eq.(26), i.e., the correction of the dispersion relation $b_p$ does not affect the rate of thermalization. On the other
hand, in the regime $b_0 \sim \Delta\beta$, the correction $b_p$ leads to a deceleration of the rate of thermalization and condensation.

## 2.4 Kinetic equation with strong disorder

In this section we discuss wave turbulence in MMFs in the presence of a 'strong' disorder, that is a random coupling among non-degenerate modes of the fiber. The regime of 'strong' disorder then differs from the previously analyzed regimes with polarization random coupling (section 2.3), which is referred to as the 'weak' disorder regime. A key distinction between the regimes of strong and weak disorder is that strong disorder prevents the conservation of beam energy during its propagation in the MMF. On the basis of the wave turbulence theory, here we derive a kinetic equation that accounts for the presence of a time-dependent random mode coupling [104]. We start the analysis by considering the NLS equation with a disordered potential that conserves the total power of the light beam, i.e., we do not consider losses. The theory describes the antagonist impacts of nonlinearity and disorder: While strong disorder enforces a relaxation to the homogeneous equilibrium distribution of the modal components ('particle' equipartition, $n_p^{eq} = \text{const}$), the nonlinear process of thermalization favours the emergence of a condensed state exhibiting a macroscopic population of the fundamental mode ($n_0 \gg n_p$ for $p \neq 0$), and an energy equipartition among the other modes. We will see that, despite a dominant strength of disorder, the system can exhibit a process of non-equilibrium condensation in the initial evolution stage, while the system eventually relaxes to the homogeneous equilibrium distribution dictated by strong disorder. The theory is confirmed by numerical simulations of the NLS equation, which are found in quantitative agreement with the simulations of the derived kinetic equation [104].

### 2.4.1 Model equation

In order to render this section self-contained, we start the presentation of the model from the original NLS equation with a random potential:

$$i\partial_z\psi = -\alpha\nabla^2\psi + V(\mathbf{r})\psi - \gamma|\psi|^2\psi + \delta V(\mathbf{r}, z)\psi. \quad (34)$$

The potential $V(\mathbf{r})$ [$\mathbf{r} = (x, y)$] is the ideal transverse index profile of the MMF, while $\delta V(\mathbf{r}, \mathbf{r})$ is the 'time'-dependent random perturbation of the potential (refractive index fluctuations), which is assumed to be real-valued with zero mean, $\langle \delta V \rangle = 0$. We recall that $\alpha = 1/(2k_0n_{co})$ and $\gamma$ denotes the nonlinear coefficient. We expand the field $\psi(z, \mathbf{r}) = \sum_p a_p(z)u_p(\mathbf{r})$ on the basis of the M real-valued eigenmodes $u_p(\mathbf{r})$ (solution of $\beta_p u_p = -\alpha\nabla^2 u_p + V(\mathbf{r})u_p$) of the unperturbed waveguide. The mode amplitudes $a_p(z)$ then satisfy

$$i\partial_z a_p = \beta_p a_p - \gamma \sum_{l,m,n} S_{plmn} a_l a_m a_n^* + \sum_l C_{pl}(z) a_p, \quad (35)$$

where the random mode coupling matrix reads

$$C_{pl}(z) = \int u_p(\mathbf{r}) \delta V(z, \mathbf{r}) u_l(\mathbf{r}) d\mathbf{r}. \quad (36)$$
The stochastic NLS Eq.(34) and the modal NLS Eq.(35) are fully equivalent. They both conserve the total power $N = \int |\psi|^2 d\mathbf{r} = \sum_p |a_p|^2$, while the random potential $\delta V(\mathbf{r}, z)$ in Eq.(34) (or the random matrix $C(z)$ in Eq.(35)), prevents the conservation of the energy (Hamiltonian).

### 2.4.2 Kinetic equation

Let us consider the situation where the random potential is a weak perturbation with respect to linear propagation effects ($\delta V \ll V$), i.e. $L_{\text{lin}} = 1/\beta_0 \ll L_{\text{dis}} = 1/\sigma_\beta$ and $L_{\text{lin}} \ll l_\beta$, where we recall that $\sigma_\beta^2$ denotes the variance of the fluctuations of the random potential (i.e., 'strength' of disorder) and $l_\beta$ denotes the corresponding correlation length. Note that this is the usual case in an optical waveguide configuration, e.g., in MMFs. Furthermore, we assume that the disorder dominates over nonlinear effects $L_{\text{dis}} \ll L_{\text{nl}}$). Note that this assumption differs from that required in the theory developed for weak disorder (see Eq.(9)), since here $L_{\text{dis}}$ does not depend on the correlation length.

Combining the (discrete) wave turbulence theory [125, 126, 127, 128, 41, 129] with tools developed for the asymptotic analysis of randomly driven linear differential equations (time-dependent disorder) [145], we derive in the Appendix C the kinetic equation governing the evolution of the averaged modal components $n_p(z) = \langle |a_p(z)|^2 \rangle$:

$$ \partial_z n_p = \sum_{l \neq p} \Gamma_{pl}^{\text{OD}} (n_l - n_p) + 8\gamma^2 \sum_{l,m,n} \frac{\delta_{plmn}^K S_{plmn}^2}{G_{plmn}^{\text{D}}} R_{plmn}[\mathbf{n}], \quad (37) $$

where

$$ R_{plmn}[\mathbf{n}](z) = n_l(z)n_m(z)n_p(z) + n_l(z)n_m(z)n_n(z) \\ - n_p(z)n_n(z)n_m(z) - n_p(z)n_n(z)n_l(z), $$

and the Kronecker symbol denotes a frequency resonance ($\delta_{plmn}^K = 1$ if $\Delta\beta_{plmn} = \beta_p - \beta_l - \beta_m + \beta_n = 0$, and zero otherwise). For clarity, we assume that the modes are not degenerate. We refer the reader to the Appendix C for the kinetic equation accounting for mode degeneracies.

The kinetic Eq.(37) unveils the interplay of nonlinearity and disorder: it reveals that diagonal and off-diagonal elements of the random matrix $C$ play fundamental different roles. The first term in the kinetic equation (37) originates in off-diagonal elements of $C_{pl}$ ($p \neq l$):

$$ \Gamma_{pl}^{\text{OD}} = 2 \int_{0}^{\infty} \langle C_{pl}(0)C_{pl}(z) \rangle \cos ((\beta_p - \beta_l)z) dz. \quad (38) $$

This equation describes an irreversible relaxation toward the homogeneous distribution, characterized by an equipartition of power ('particles') among the modes, $n_p^{\text{eq}} = N/M = \text{const}$. This process occurs over the typical propagation length

$$ \mathcal{L}_{\text{kin}}^{\text{eq}} \approx 1/\overline{\Gamma_{jl}^{\text{OD}}}. \quad (39) $$
This is the well-known evolution of a system ruled by strong random mode coupling. To be more precise, we note that $\mathcal{L}_{\text{kin}}^{\text{eq}} = 1/\lambda_2^{\text{OD}}$ where $-\lambda_2^{\text{OD}}$ is the second eigenvalue of the matrix $\tilde{\Gamma}^{\text{OD}}$, with $\tilde{\Gamma}_{jl}^{\text{OD}} = \Gamma_{jl}^{\text{OD}}$ for $j \neq l$ and $\tilde{\Gamma}_{jj}^{\text{OD}} = -\sum_{l \neq j} \Gamma_{jl}^{\text{OD}}$, the first eigenvalue being $\lambda_1^{\text{OD}} = 0$.

We now show that this relaxation process mediated by strong disorder does not necessarily inhibit the nonlinear processes of thermalization and condensation. This becomes apparent through the analysis of the collision term in the kinetic equation (37), which exclusively involves the diagonal components $C_{pp}$:

$$ \Gamma_{pl}^{\text{D}} = \int_{0}^{\infty} \langle C_{pp}(0)C_{ll}(z) \rangle + \langle C_{ll}(0)C_{pp}(z) \rangle dz. \quad (40) $$

The matrix $\Gamma^D$ contributes to the tensor involved in the collision term (see the Appendix C):

$$ G_{plmn}^{\text{D}} = \Gamma_{ll}^{\text{D}} + \Gamma_{mm}^{\text{D}} + \Gamma_{nn}^{\text{D}} + \Gamma_{pp}^{\text{D}} + 2\Gamma_{lm}^{\text{D}} - 2\Gamma_{ln}^{\text{D}} - 2\Gamma_{lp}^{\text{D}} - 2\Gamma_{mn}^{\text{D}} - 2\Gamma_{mp}^{\text{D}} + 2\Gamma_{np}^{\text{D}}. $$

To discuss the role of the collision term, let us forget for a while the first term in the kinetic equation (37). The collision term conserves $N = \sum_p n_p(z)$, $E = \sum_p \beta_p n_p(z)$, and exhibits a H-theorem of entropy growth $\partial_z S \ge 0$, with $S(z) = \sum_p \log[n_p(z)]$, see Appendix C. Hence, the collision term describes a process of thermalization to the RJ distribution $n_p^{\text{RJ}} = T/(\beta_p - \mu)$, exhibiting energy equipartition among the modes. This occurs over a typical propagation length

$$ \mathcal{L}_{\text{kin}}^{\text{RJ}} \approx L_{nl}^2 \overline{G_{plmn}^{\text{D}} / S_{plmn}^2}. \quad (41) $$

As discussed in the Appendix B, for an energy smaller than the critical value $E \le E_{\text{crit}}$, the RJ distribution $n_p^{\text{RJ}}$ exhibits a phase transition to a condensed state [12]. The condensate amplitude $n_0$ then constitutes the natural parameter that distinguishes the two antagonist regimes:

(i) For $\mathcal{L}_{\text{kin}}^{\text{eq}} \ll \mathcal{L}_{\text{kin}}^{\text{RJ}}$, the disorder dominates and $n_0(z) \to n_p^{\text{eq}} = N/M = \text{const}$ for $p = 0, 1, \dots, M-1$;

(ii) For $\mathcal{L}_{\text{kin}}^{\text{eq}} \gg \mathcal{L}_{\text{kin}}^{\text{RJ}}$, the dynamics is dominated by RJ thermalization, and condensation leads to a macroscopic population of the fundamental mode $n_0(z) \to n_0^{\text{RJ}} \gg n_p^{\text{RJ}}$ for $p = 1, \dots, M-1$.

### 2.4.3 Numerical simulations

We have performed extensive numerical simulations to test the validity of our theory. We have considered the concrete example of a parabolic trapping potential of the form $V(\mathbf{r}) = q_x x^2 + q_y y^2$, with the fundamental mode eigenvalue $\beta_0 = \sqrt{\alpha}(\sqrt{q_x} + \sqrt{q_y})$. Let us consider a general model of disorder with a random potential of the form $\delta V(\mathbf{r}, z) = \mu(z)g(\mathbf{r})$, where $\mu(z)$ is a real-valued stochastic function with zero mean and $\langle \mu(0)\mu(z) \rangle = \sigma_{\beta}^2 \exp(-|z|/l_{\beta})$. In order to remove mode degeneracies, we consider in the simulations an elliptical parabolic potential ($q_x \neq q_y$). To compute the matrices $\Gamma^{\text{OD}}$ and $\Gamma^{\text{D}}$ in analytical form, we consider the random mode coupling $g(x, y) = \cos(b_x x/x_0) \cos(b_y y/y_0)$, where $(x_0, y_0)$ denote the radii of the fundamental elliptical mode, see Appendix C.
<figure><img src="image_9.png" /><figcaption>Figure 8: Dynamics dominated by strong disorder $\mathcal{L}_{\text{kin}}^{\text{RJ}} \gg \mathcal{L}_{\text{kin}}^{\text{eq}}$. The system irreversibly relaxes toward the equilibrium $n_p^{\text{eq}}$. Evolutions of the fundamental mode $n_0(z)$ (a), and $n_2(z)$ (b), the energy $E(z)/(N\beta_0)$ (c), obtained from the numerical simulation of the NLS Eq.(35): 64 realizations are reported with colored lines; the bold white line is the corresponding empirical average; the dashed black line is the prediction of the KE (37). (d) Modal distribution $n_p$ in the initial condition ($z=0$, blue) and at $z\beta_0 = 2 \times 10^6$ for the NLS simulation (red), and the kinetic equation (black). Parameters: $L_{\text{dis}}/L_{\text{lin}} = 7$, $L_{\text{dis}}/L_{\text{nl}} = 4.1 \times 10^{-4}$, $l_\beta\beta_0 = 42$. Source: From Ref.[104].</figcaption></figure>
<figure><img src="image_10.png" /></figure>

Figure 9: Thermalization precedes equilibrium relaxation: Same panels as in Fig. 8, but in the regime $\mathcal{L}_{\text{kin}}^{\text{RJ}} \lesssim \mathcal{L}_{\text{kin}}^{\text{eq}}$. The system exhibits an incipient process of RJ thermalization and nonequilibrium condensation characterized by a growth of the condensate amplitude $n_0(z)$ for $z\beta_0 \lesssim 2 \times 10^5$. Disorder subsequently prevails, which induces a decay of $n_0(z)$ (and eventually brings the system to equilibrium $n_p^{\text{eq}}$). Parameters: $L_{\text{dis}}/L_{\text{lin}} = 7$, $L_{\text{dis}}/L_{\text{nl}} = 0.033$, $l_{\beta}\beta_0 = 167$. Source: From Ref. [104].
<figure><img src="image_11.png" /><figcaption>Figure 10: Thermalization prevails over equilibrium relaxation: Same panels as in Fig. 8, but in the regime $\mathcal{L}_{\text{kin}}^{\text{RJ}} \ll \mathcal{L}_{\text{kin}}^{\text{eq}}$. The system exhibits a process of RJ thermalization and condensation characterized by a significant growth of the condensate amplitude $n_0(z)$ to the value predicted by the RJ distribution, $n_0^{\text{RJ}}/N \simeq 0.55$ (horizontal dashed-dotted black line) (a). At variance with Figs. 8-9, the energy $E(z)$ is almost constant (c). The modes approach the RJ distribution $n_p^{\text{RJ}}$ (green) (d). Parameters: $L_{\text{dis}}/L_{\text{lin}} = 8.4$, $L_{\text{dis}}/L_{\text{nl}} = 0.04$, $l_{\beta}\beta_0 = 4 \times 10^3$. Source: From Ref.[104].</figcaption></figure>
According to the theory, the two terms in the kinetic equation (37) are an-
tagonists and compete against each other. If $\mathcal{L}_{\text{kin}}^{\text{eq}} \ll \mathcal{L}_{\text{kin}}^{\text{RJ}}$, disorder prevails and
the system relaxes to the expected equilibrium $n_p^{\text{eq}}$ =const. This is illustrated in
Fig. 8, which reports the results of the numerical integration of the NLS Eq.(35)
for 64 realizations ($\mathcal{L}_{\text{kin}}^{\text{RJ}}/\mathcal{L}_{\text{kin}}^{\text{eq}} \simeq 400$). The corresponding average over such
realizations (bold white line) is in agreement with the simulation of the kinetic
equation (37) (dashed black line) starting from the same initial condition. Here
and thereafter, the quantitative agreement between NLS and kinetic equation
simulations is obtained without any adjustable parameter.

Unexpectedly, however, a nonequilibrium process of condensation and ther-
malization can be observed in the initial stage of propagation when $\mathcal{L}_{\text{kin}}^{\text{RJ}} \lesssim \mathcal{L}_{\text{kin}}^{\text{eq}}$
(see Fig. 9 for $\mathcal{L}_{\text{kin}}^{\text{RJ}}/\mathcal{L}_{\text{kin}}^{\text{eq}} \simeq 0.07$), while asymptotically the system still relaxes
to the homogeneous equilibrium state $n_p^{\text{eq}}$. The nonequilibrium property of
condensation is reflected by the fact that the energy $E(z) = \sum_p \beta_p n_p(z)$ is not
conserved during the evolution, see Fig. 9(c).

We stress that the condensation processes may occur very efficiently by in-
creasing the correlation length $l_\beta$, in such a way that $\mathcal{L}_{\text{kin}}^{\text{RJ}} \ll \mathcal{L}_{\text{kin}}^{\text{eq}}$, see Fig. 10 for
$\mathcal{L}_{\text{kin}}^{\text{RJ}}/\mathcal{L}_{\text{kin}}^{\text{eq}} \simeq 0.003$. In this regime, the energy is almost conserved $E \simeq \text{const}$ and
RJ thermalization occurs almost completely, as confirmed by the modal popu-
lations that approach the RJ distribution $n_p^{\text{RJ}}$ (Fig. 10(d)), and the condensate
approaches the RJ prediction $n_0^{\text{RJ}}/N \simeq 0.55$, see Fig. 10(a). Note that a similar
effect of pre-thermalization and condensation has been identified recently in
nonlinear disordered Floquet systems [146].

It is interesting to remark that, although random mode coupling does not
affect the RJ equilibrium, it significantly accelerates the rate of thermalization
to this equilibrium. Indeed, the ratio of thermalization in the presence and
the absence of disorder scales as ~ $\bar{G}^D/\beta_0 \ll 1$. This effect of acceleration
of thermalization is very similar to that pointed out in the previous section
through the analysis of weak disorder, see Sec. 2.3.6.

2.5 *Conclusions of the wave turbulence approach*

On the basis of the wave turbulence theory, we have derived kinetic equa-
tions describing the nonequilibrium evolution of random waves in a regime
where disorder dominates nonlinear effects ($L_d \ll L_{nl}$). The theory revealed
that the presence of a conservative weak disorder introduces an effective dissi-
pation in the system, whose resonance broadening prevents the conservation
of energy, which inhibits the effect of condensation within the usual continuous
wave turbulence approach. However, usual experiments of beam-cleaning are
not described by the continuous wave turbulence theory, but instead by a dis-
crete wave turbulence approach. In this discrete turbulence regime only exact
resonances contribute to the kinetic equation, which is no longer sensitive to
the effect of dissipation-induced resonance broadening. Accordingly, the dis-
crete kinetic equation conserves the energy, which re-establishes the process of
wave condensation. The main result is that the effective dissipation induced by
disorder modifies the regularization of such discrete resonances, which leads
to an acceleration of the rate of thermalization and condensation.
We have considered different models of weak disorder in MMFs. The theory reveals that when all modes experience the same (mode-correlated) noise, dissipation induced by disorder vanishes, and the system no longer exhibits a fast process of condensation. However, even a relative small decorrelation among the noise experienced by the modes is sufficient to re-establish a disorder-induced acceleration of condensation. The simulations are in quantitative agreement with the theory without using adjustable parameters, in a regime where disorder dominates nonlinear effects ($L_d \ll L_{nl}$). However, the impact of weak disorder due to polarization fluctuations in beam cleaning experiments is expected to be of the same order as that of nonlinear effects [118]. Accordingly, the impact of a moderate disorder ($L_d \gtrsim L_{nl}$) was considered in Ref.[25]. In this regime, the phase-correlations among the modes are not completely suppressed, and the system enters a mixed coherent-incoherent regime. Accordingly, the modal components can exhibit an oscillatory behavior that slows down the thermalization process, although the optical beam still exhibits a fast condensation process that is consistent with the experimental results of RJ thermalization. In addition, the analysis of polarization effects reported in [25] revealed that optical beam cleaning is responsible for an effective partial repolarization of the central part of the beam. This property was observed experimentally in Ref.[147] and it can be interpreted as a consequence of the condensation-induced macroscopic population of the fundamental mode of the MMF. Then considering the dominant contribution of weak disorder originating in polarization fluctuations, the nonequilibrium kinetic theory and simulations both provide a qualitative understanding of the effect of optical beam self-cleaning, in particular when a large number of modes are excited into the MMF, as in Ref.[12]. On the other hand, the discrete nature of wave turbulence in MMFs is responsible for an effective freezing of thermalization and condensation when only a small number of modes is excited. Such an effective freezing of condensation may also explain why beam cleaning has not been clearly observed in step-index MMFs [25].

It is important to recall that, at variance with the experiments of beam cleaning where sub-nanosecond pulsed are injected into the MMF, in this section we considered a purely spatial model of light propagation in the fiber. This is justified by the fact that in the nanosecond regime temporal dispersion effects play a negligible role at leading order. Then although the theory can describe a mechanism of beam cleaning condensation, it cannot provide a complete description of the experimental results.

In the second part of this section we have derived a wave turbulence kinetic equation in the presence of strong disorder. The theory revealed that RJ thermalization and condensation can even take place when the dynamics is dominated by strong disorder ($L_{dis}/L_{nl} \ll 1$), provided that $\mathcal{L}_{\text{kin}}^{\text{RJ}} \ll \mathcal{L}_{\text{kin}}^{\text{eq}}$. Note that, at variance with weak disorder, strong disorder introduces very large fluctuations of the modal components $n_p$ (see NLS simulations in Figs. 8-10). The comparison with the theory then require an average over the realizations, which is found in quantitative agreement with the simulations of the kinetic Eq.(37), without using adjustable parameters.
Experiments carried out using relatively short pieces of MMFs (~12m) have been reported in Ref.[104] to test the theory at a purely qualitative level. Strong mode coupling was obtained by applying a stress to the MMF with clamps [148]. By adjusting the applied stress, the strength of mode coupling could be tuned. Strong random mode coupling leads to an increase of the energy (Hamiltonian) – the larger the strength of applied stress, the larger the increase of energy. The experiments show that in the presence of significant applied stress, RJ thermalization and condensation are inhibited. By decreasing the amount of applied stress, the experiments confirm that RJ condensation may still occur, even in the presence of strong disorder [104].

It is important to note that an effect of beam cleaning of different nature than RJ condensation has been recently observed experimentally in MMFs in the linear regime of propagation, over very large propagation distances of the order of the kilometer, see Ref.[149]. With such large propagation lengths, strong random mode coupling in MMFs can no longer be neglected. Taking into account for the losses due to a linear coupling between the highest order guided modes and the radiative modes, a non-homogeneous solution predicting a larger population of the fundamental fiber mode was obtained [150] (also see [151]). This can provide a natural explanation for the remarkable effect of beam cleaning observed experimentally in the linear regime in Ref.[149].

The linear propagation of scalar waves in randomly perturbed open waveguides was also addressed in Ref.[152]. From first principles and by using a multiscale analysis, it was shown that the coupling between guided and radiative modes involves a perturbed form of 'particle' equipartition: the mean powers of the guided modes become proportional to the first eigenvector of an effective scattering matrix that is a refined version of (38), and that is determined by the second-order statistics of the waveguide fluctuations. When the coupling from guided to radiative modes is negligible compared to the coupling between guided modes, this eigenvector is uniform and the standard homogeneous equipartition result is recovered. But in general the eigenvector is not uniform, and it can exhibit a concentration close to the lowest eigenmode, for which the coupling towards the radiative modes is the smallest.

### 3 Experimental evidences

In this Section we present a series of recent experimental studies, aiming at validating the predictions of the thermodynamic or statistical mechanics approach. In these works, pulsed lasers with peak powers of a few tens of kW with pulse durations ranging between hundreds of femtoseconds to nanoseconds were used. In these conditions, it has been verified that thermodynamic equilibrium can be reached for beams propagating over the relatively short distances of a few meters of standard GRIN MMF. With such power levels and propagation distances, one fulfills the theoretical hypotheses of conservation of the optical beam power (or 'number of particles') and kinetic energy upon propagation in the MMF, given that both linear and nonlinear optical losses can be neglected [153, 154], as discussed in section 2. However, working with short
pulses may lead to additional nonlinear effects, that modify their temporal and spectral properties. The main effects here are self-phase modulation, Raman scattering, and geometric parametric instability, which produce a progressive beam spectral broadening upon propagation. In our statistical description of BSC, we have considered purely monochromatic waves, disregarding temporal aspects. To alleviate this, experimental measurements were typically carried out either by integrating pulses over time or, as an alternative, by inserting at the fiber output a narrow bandpass filter centered at the laser wavelength. In this regard, it is important to mention that Mangini et al. showed that the RJ distribution at thermal equilibrium may be reached, for input pulse durations ranging from hundredths of ps down to hundredths of fs (see Fig. 11).

<figure><img src="image_12.png" /><figcaption>Figure 11: Illustration of the mode power distribution when sorting the modes by their radial (ℓ) and azimuthal (m) indices, for input pulse duration of 174 fs (a), 1 ps (b), 7.6 ps (c), and 435 ps (d). The black bars are the experimental values, while the red lines are the theoretical RJ distributions. Image from [108].</figcaption></figure>

Direct observations of the attainment of the RJ equilibrium were made possible by the development of holographic MD tools. A notable example is that of MD systems based on phase-only spatial light modulators, which allow for fully decomposing a beam in its modal components, i.e., both the mode amplitude and phase can be retrieved [155]. Other approaches exploited, instead, interferometric setups, e.g., in Ref. [107]. Generally speaking, holographic methods were quite successful in experimental studies of the statistical mechanics theory.

Within this experimental framework, many investigations of the BSC effect were reported in the literature [107, 108, 109, 110, 111, 156, 157, 158]. Still, a few papers studied the thermalization process of spatially structured light (see Ref.
[141] for a comparison between wave thermalization and BSC). In particular, some experiments have shown that, when increasing the beam power well above the self-cleaning threshold, a thermalized beam may still yield a speckled intensity pattern [156]. In addition, the presence (and conservation) of the orbital angular momentum of light may lead to thermalized beams without a bell shape [26]. In this regard, in Fig. 11 one may note that each couple of bars corresponding to modes with the same value of $l$, and opposite signs of $m$ have the same height. Furthermore, as an ultimate case, whenever a beam is associated with a negative temperature, the occupancy of the fundamental mode tends to vanish, thus making it virtually impossible to achieve a Gaussian-like intensity profile [141].

When it comes to BSC, it may seem counterintuitive that as the beam cleans up, its associated entropy increases, i.e., the disorder grows larger. This has led to speculations that the whole statistical mechanics approach to BSC fails because of a breach of the second law of thermodynamics [159]. Nevertheless, the growth of entropy during BSC was observed in several experiments. Baudin et al. have determined the equilibrium entropy in the classical RJ condensation in a GRIN MMF, proving that the equilibrium entropy decreases with the kinetic energy [12].

Still, verifying that BSC is consistent with the second principle of thermodynamics means verifying that the entropy reaches a maximum value when thermal equilibrium is attained. This poses the challenges of experimentally measuring the entropy of out-of-equilibrium states.

In the following, we will go through recent experimental results that have achieved such a goal [109, 110]. Then, we will present the cases of wave thermalization in three significant cases, i.e., at high, near zero, and negative temperatures. Finally, we will shortly review the recent discovery of a modal phase-locking mechanism which accompanies BSC.

## 3.1 The growth of optical entropy

A remarkably simple, and yet effective way of writing the Boltzmann entropy was proposed by Wu et al. in [13]. This reads

$$S = \sum_{p} \ln n_{p}, \qquad (42)$$

where the sum runs over all guided modes. Strictly speaking, such a formula only holds at thermal equilibrium. However, it is interesting to note that the same expression of the entropy also holds out of equilibrium, where the modal populations depend on the propagation length, as discussed through the H-theorem of entropy growth in the previous section 2, see Eqs.(28). Accordingly, the nonequilibrium entropy grows upon beam propagation, and its growth saturates as the field approaches thermal equilibrium, where the entropy reaches its equilibrium value. In this way, the entropy given in Eq.(42) (with $n_p$ given by the RJ distribution) provides the upper bound of the out-of-equilibrium entropy. As discussed in section 2, this occurs under the hypothesis that the fiber modes have Gaussian statistics, which may result from
the presence of random linear mode coupling (i.e., fluctuations of the index of refraction).

According to (42), to experimentally verify the entropy growth one needs to track the evolution of the occupation of the modes, i.e., the values of the {$n_p$}, over "time". As discussed above, the role of time in nonlinear wave thermalization is played by the propagation distance $z$. However, the BSC effect is typically observed in standard GRIN MMFs by increasing the beam input power at a given fiber length, with fixed injection conditions. This is mostly because varying the power is a trivial experimental operation, whereas cut-back experiments where the fiber is progressively shortened may be difficult and tedious (especially when short fiber spans are used). In addition, one needs to realign the output fiber tip at each cut, which may be a tough process when operating with a holographic MD setup. Finally, it has to be noted that each cutback experiment requires a brand-new fiber. Interestingly, a specular problem is faced when running simulations: it is quite simple to add steps of integration of the nonlinear wave equation; whereas quasi-continuously varying the beam power requires running several simulations, thus being quite time-consuming.

In standard experiments, increasing the beam power has the beneficial effect of facilitating its thermalization. This is because higher powers provide enhanced nonlinear effects, i.e., more efficient power and energy transfers among the modes. Indeed, in the thermodynamic-analogy picture, one may see non-linearity as the mechanism leading to interactions among the modes via photon exchanges. Now, from (42) it is evident that the entropy logarithmically scales with power. Indeed, if one doubles the power of all the modes, the entropy is enhanced by a $M \ln 2$ factor, being $M$ the number of modes. Therefore, it is convenient to separate the explicit and the implicit contributions of power to the entropy. This can be easily done by rewriting the Boltzmann entropy as

$$S = M \ln N + \sum_{p=1}^{M} \ln |f_p|^2, \qquad (43)$$

where $|f_p|^2 = n_p/N$ is the power fraction of the p-th mode. In fact, in (43) the first term contains the logarithmic increases of $S$ with $N$, while the second term only implicitly depends on $N$. Within this expression, it appears evident that to prove the entropy growth during thermalization one can just compute the second term in (43), which is referred to as configuration entropy [109] (or entropy per particle in [110]),

$$\tilde{S} = \sum_{p=1}^{M} \ln |f_p|^2. \qquad (44)$$

Indeed, no matter whether wave thermalization is achieved by increasing the power at a fixed fiber length, or at a fixed input power and then increasing the length as it can be verified by cut-back experiments, the maximization of $\tilde{S}$ directly proves the maximization of $S$.
A recent work has reported both types of experiments [109]. The main results are illustrated in Fig. 12. The value of $\tilde{S}$ were calculated by eq.(44) starting from the values of $n_p$ that were experimentally measured with the holographic MD method of Ref. [155]. In the first set of experiments, whose results are shown in Fig. 12a, a trend of growing $\tilde{S}$ at a fixed length was observed, when increasing the input power above the self-cleaning power threshold. For a short fiber length (2 m), a marked increase of $\tilde{S}$ is observed when the power grows up to 5 kW (cfr. blue curve in Fig. 12a). At higher powers, $\tilde{S}$ shows a saturation trend. This is expected when the beam is near thermal equilibrium, since the entropy has to reach its maximum. Interestingly, Fig. 12a also shows that at a longer fiber length (12 m), $\tilde{S}$ experiences, instead, a relatively slow growth (cfr. red curve in Fig. 12a). This is because the output beam remains relatively close to its state of equilibrium at all input power values.

The results of the second type of experiments (cutback) are shown in Fig. 12b: as can be seen, here $\tilde{S}$ exhibits a general trend of growing larger with distance. Consistently with the results of Fig. 12a, at the highest peak power of 14.9 kW, $\tilde{S}$ remains essentially a constant for lengths longer than 4 meters. This can be explained by considering that at these power levels thermalization is achieved at relatively shorter distances. On the other hand, an overall rapid growth of $\tilde{S}$ over the first 6 meters is observed at lower powers (see, for instance, the light blue curve in Fig. 12b). Whenever thermal equilibrium is reached, the entropy growth ends with a plateau.

So far, we have considered the growth of entropy of a single beam, i.e., its transition through non-equilibrium states to attain a final state thermal equilibrium. In a loose sense, this process can be interpreted as a consequence of the second law of thermodynamics, which has a probabilistic origin in the statistical mechanics framework. As such, thermodynamic parameters such as the temperature and the chemical potential merely acquire a statistical meaning. It is, indeed, when heat exchanges between two systems are allowed that one may provide $S, T$, and $\mu$ with a physical connotation, without recurring to probabilistic arguments. Indeed, the entropy growth can be seen as the driver of the unidirectional nature of the heat flow from hot to cold. Hence, in accordance with the zeroth principle of thermodynamics, one may demonstrate that $S, T$, and $\mu$ are true physical forces by observing an irreversible energy transfer between two systems with different temperatures via irreversible heat exchanges. This result was experimentally achieved in Ref. [110], which will be briefly reviewed in the following subsection.

### 3.1.1 Two-beam thermalization and optical calorimetry

Thermodynamics involves the description of the properties of physical systems at thermal equilibrium via macroscopic parameters such as temperature [83]. On the other hand, calorimetry deals with thermodynamic processes under non-equilibrium conditions, e.g., via the exposition of a given object to an external heat source. The easiest calorimetry experiment that one can imagine involves two bodies with different temperatures, that are put in contact inside
<figure><img src="image_13.png" /><figcaption>Figure 12: a) Evolution of configuration entropy vs. input peak power during BSC for a fiber length of 2 (blue symbols) and 12 m (red symbols). b) Increase of the configuration entropy vs. fiber length which accompanies BSC, for different input beam powers (the legend shows the average and the peak power of the optical pulses, respectively). Symbols are experimental data, while the dashed curves are mere guides for the eye. The experiments were carried out with a standard 50/125 GRIN MMF, coiled around a 15 cm diameter plastic drum. Figure from [109].</figcaption></figure>

a calorimeter. In the latter, heat exchanges are allowed between the bodies, but there is no heat dissipation towards the outside world. In other words, the system composed of the two objects is adiabatic. Within this setup, the heat uni-directionally flows from the higher-temperature body towards the other body. An equilibrium is then established whenever the two bodies reach the same temperature, and heat exchanges no longer take place.
Such an experiment was recently reproduced within the multimode pho-
tonics platform. Two laser beams, representing the two bodies, were simul-
taneously injected into an MMF. The two beams are distinguished by their
polarization state. The first beam (A), having a temperature $T_A$ and chemical
potential $\mu_A$, is horizontally polarized; whereas the second beam (B), which is
hotter ($T_B > T_A$, $\mu_B > \mu_A$), has a vertical state of polarization. The initial state
of the system associated with this experimental setup can be depicted as two
gases that occupy two halves of the total available mode volume (see Fig. 13a).
On each half of the box, the color of the beams indicates their temperature.
Specifically, the colder (hotter) gas is depicted in blue (red). This use of colors,
however, is only meant to help the visualization of the phenomenon. Indeed,
as well known, photons are indistinguishable particles.

It has to be noted that, in analogy with the previously discussed case of
the entropy, it is possible to formally associate a beam with both an optical
temperature and chemical potential even when it is not yet at equilibrium.
Indeed, one may compute the values of T and μ starting from the knowledge of
the mode occupancies, even when these do not follow the RJ law, as depicted at
the bottom of Fig. 13a. This is, strictly speaking, improper since the temperature
and chemical potential do not possess a rigorous physical meaning whenever
one is out of equilibrium. In this sense, one may refer to T and μ as the values of
temperature and chemical potential that a beam would eventually reach, once
it is at thermal equilibrium. In the case of two beams considered here, T and μ
of each beam are those achieved in the absence of the other beam, i.e., with a
single beam injection into the fiber.

The initial state in Fig. 13a is defined by $T^H = T_A$, $\mu^H = \mu_A$, $T^V = T_B$, and $\mu^V = \mu_B$, where the superscript letters 'H' and 'V' are used to identify the parameters of the left and right sides of the box, respectively. Once the system evolves, the two photon gases interact, exchanging both energy and particles (cfr. Fig. 13b), until progressively reaching the same equilibrium values of temperature ($T_{eq}$) and chemical potential ($\mu_{eq}$), as depicted in Fig. 13c. The exchange of particles is due to the lack of polarization conservation during single beam propagation. Indeed, it was shown in Ref. [157] that, within standard experimental conditions, an optical beam undergoes full depolarization.

The equilibrium values of temperature ($T_{eq}$) and chemical potential ($\mu_{eq}$) are first established at the center of the box, where the initial interactions between the two gases occur (refer to Fig. 13b). For the sake of simplicity, we indicate the equilibrium region by depicting violet particles in Figs. 13b-c. Thermal equilibrium is then gradually achieved at the periphery of the box (see Fig. 13c). Note, however, that the boxes in Fig. 13 indicate a mode volume. Thus, defining a "local" temperature is formally improper. During the nonequilibrium transient state, each side of the box has a mode distribution that gradually approaches the RJ law (as depicted in the histograms at the bottom of Figs. 13b,c). The color map indicates that before equilibrium is reached, the temperatures on the two sides of the box differ. The left side of the box still contains the colder gas, resulting in a higher temperature for the horizontally polarized photon gas when compared with the vertically polarized photons, as the for-
mer still contains more energetic (red) particles (i.e., $T_A < T^H \le T^V < T_B$, and $\mu_A > \mu^H \ge \mu^V > \mu_B$). The temperature difference between the left and right sides of the box diminishes as the input power increases, eventually disappearing when thermal equilibrium is reached (i.e., $T^H = T^V = T_{eq}$ and $\mu^H = \mu^V = \mu_{eq}$).

This calorimetry picture was confirmed by experiments that rely on the same holographic MD technique used in the case of a single beam (see Ref. [110] for details). In particular, the growth of the optical entropy was experimentally verified, in pure analogy with the results presented in Sec. 3.1. Here, however, the relative power of the two beams (i.e., the ratio of the mass of the two photon gases) provided an additional degree of freedom. Specifically, the mass ratio impacts the final equilibrium temperature. Moreover, the latter can be varied by tuning the temperature of only one of the two beams. This result, which was experimentally demonstrated, is of utmost importance for applications based on BSC. Indeed, the success of optical calorimetry experiments showed that it is possible to modify the optical temperature, i.e., the quality or brightness of a laser beam, by exploiting the nonlinear interaction with another beam in an MMF. In this sense, these results pave the way for developing a new generation of photonic tools for all-optically controlling the spatial profile of intense laser beams. The fact that the temperature of a beam is associated with its spatial quality was further highlighted in Ref. [156], whose content will be shortly reviewed below.

## 3.2 High-temperature wave thermalization

In Sec. 3.1 and, specifically, in Fig. 12, we have discussed the establishment of thermal equilibrium when increasing either the input power $P$ or the fiber length. However, as we have previously mentioned, these two methods are not equivalent. Indeed, on the one hand, the propagation distance plays the role of time. Therefore, the number of particles and the kinetic energy are conserved in cutback experiments: thus, one may properly associate the beam with a temperature $T$ and a chemical potential $\mu$. On the other hand, when achieving wave thermalization by increasing the input power, the number of particles and the kinetic energy are no longer conserved. Thus, the values of $T$ and $\mu$ continuously vary with $P$. Specifically, increasing $P$ leads to an increase of $T$, i.e., to a larger occupancy of the HOMs. As a result, too high input powers cannot provide BSC, even within the weak nonlinearity regime.

This concept was experimentally demonstrated in Ref. [156], by using the same holographic MD setup of the experiments discussed in Sec. 3.1. The main result is illustrated in Fig. 14. Specifically, the insets show the output transverse intensity profiles obtained for increasing values of the beam power. When the input peak power grows from 1 kW up to about 13 kW, the BSC effect takes place: the beam transforms from speckles to a bell shape. However, as the input power is further increased, a clear "self-spoiling" of the output beam occurs.

A quantitative analysis of the intensity patterns ($I$) shown in Fig. 14 was carried out by computing their correlation with the intensity of the fundamental
<figure><img src="image_14.png" /><figcaption>Figure 13: Sketch of the thermalization process of two photon gases. (a) Initial state: two orthogonally polarized beams are prepared with different temperatures, chemical potentials, and mode distributions. (b) Nonequilibrium transient state: the difference in temperature and chemical potential among the left and right sides of the box progressively quenches, until their associated mode occupancy progressively approaches the RJ distribution, as a consequence of nonlinear beam propagation, i.e., of heat exchange between the two gases. (c) Equilibrium state: the two gases reach the same values of temperature and chemical potential, and their mode distribution obeys the RJ law. Figure from [110].</figcaption></figure>

mode ($I_0 \propto \exp\{-2(x^2 + y^2)/w^2\}, w = 6.33 \mu m$) as

$$Cor = \frac{\int I(x, y) I_0(x, y) dxdy}{\sqrt{\int I(x, y)^2 dxdy \int I_0(x, y)^2 dxdy}}, \quad (45)$$
The experimental variation of *Cor* with input beam power is reported in the central panel of Fig. 14. A maximum value of *Cor* ≃ 98% is reached at the occurrence of BSC. On the contrary, at either higher or lower powers, *Cor* rapidly drops. An analogous (and opposite) trend was observed for the *M*² parameter (not shown here, see Ref. [156]). In the linear regime, a value of *M*² ≃ 3 was measured. Whereas, near the BSC threshold, *M*² decreases below 2. Then, at higher powers, the beam quality worsens, and the associated *M*² grows even more than in the linear regime.

In this regard, it is important to mark the difference between the speckled
patterns at low and high powers. The former corresponds to nonequilibrium
states, whereas the latter are associated with thermalized waves at high tem-
peratures. This was confirmed by MD experiments: the distribution of the
mode occupancy at low input power only depends on the beam injection con-
ditions; to the contrary, the degradation of the beam quality at powers above
the BSC threshold, i.e., the beam self-spoiling, always leads to a mode power
distribution that obeys the RJ law.

From the application point of view, these findings establish a maximum
brightness for multimode beam propagating in nonlinear GRIN MMF, at least
in the weakly nonlinear regime.

In summary, high optical temperatures lead to the loss of spatial beam
quality, which is detrimental for most fiber-based applications. In contrast,
as discussed in Sec. 3.1.1, cooling down an optical beam increases the power
fraction of the fundamental mode. The ultimate case would be *T* = 0, i.e.,
the full wave condensation condition in which a single mode beam propagates
undisturbed in an MMF. In contrast, whenever *T* < 0 the fundamental mode is
depleted in favor of HOMs. Both the *T* near-zero and the negative temperature
cases will be the topics of the next section.

3.3 *Light thermalization and condensation by varying the beam energy*

In this section, we discuss experiments in MMFs that have been performed
by varying the kinetic energy of the launched optical beam (we remind here
that the kinetic energy refers to the linear contribution of the Hamiltonian,
see section 2.1). In this way, we present two different phenomenologies, the
RJ condensation of classical waves, and the thermalization toward negative
temperature equilibrium states.

3.3.1 Rayleigh-Jeans condensation

In spite of the physical differences that can distinguish RJ condensation
and BE condensation of quantum particles (atoms, molecules, polaritons, pho-
tons...), the underlying mathematical origin is similar, because of the common
singular behavior (vanishing denominator) of the equilibrium Bose distribution
for quantum particles, and the equilibrium Rayleigh-Jeans (RJ) distribution for
classical waves, as recently discussed in Ref.[160]. The RJ condensation of clas-
sical waves has been studied since a long time [161, 162, 57, 11], and we refer the
reader to the Appendix B for a succinct summary of the equilibrium properties
of the phase transition of wave condensation.
<figure><img src="image_15.png" /><figcaption>Figure 14: Center panel: correlation between the output beam intensity profile and the fundamental GRIN fiber mode, vs. input peak power; insets show the experimental output beam intensity for increasing values of the peak power *P* of 7.8 ps input pulses, showing that the initial beam cleaning (green background region) is followed by beam spoiling at higher input power values (violet background region). Source: From Ref. [156].</figcaption></figure>

As commented in previous sections in relation with RJ thermalization, light propagation in MMFs is essentially a conservative system, in which the power ('particle number') and the kinetic energy are conserved through propagation. In other terms, there is no thermostat in the optical experiments: By keeping constant the power ("number of particles") *N*, the energy *E* plays the role of the control parameter in the transition to condensation for this microcanonical statistical ensemble. Here, we discuss the experimental observation of the transition to condensation by varying the energy *E* of the beam launched into the MMF, while keeping constant the power *N* [12]. In order to vary the kinetic energy, the laser beam is passed through a diffuser to generate a speckle beam, which is subsequently injected in the fiber. In this way, the amount of randomness of the speckle beam can be varied by keeping a constant
power, i.e., the larger the randomness of the speckle beam, the higher the energy. Following this procedure, the experimental results are averaged over an ensemble of different realizations of speckle beams with the same power and the same energy. The importance of this aspect can be understood by recalling the fact that the RJ distribution is inherently a statistical distribution, whose verification requires, in principle, to consider an average over a statistical ensemble of speckle beams.

The experiments reported in Ref.[12] have been performed by using a conventional GRIN MMF of 12 m long that guides $M \simeq 120$ modes (fiber radius $R = 26$ µm). In the experiments, the averages over ensembles of speckle beams have been performed from a total number of $2 \times 1000$ measurements of the near-field and far-field intensity distributions recorded for different energies $E$ and a fixed power of $N = 7$ kW. The ensemble of realizations of speckle beams was partitioned within small energy intervals $\Delta E$, and then averaged in each individual energy interval, see Ref.[12] for details.

At variance with homogeneous condensation in a bulk medium (i.e., without a trapping potential $V(r) = 0$) [57, 41], the presence of the parabolic potential reestablishes condensation in the 'thermodynamic limit' in 2D. There exists a (non-vanishing) critical energy $E_{crit}^* = NV_0/2$ such that $\mu = \beta_0$, where $V_0$ denotes the depth of the parabolic trapping potential associated with the finite number of modes of the GRIN MMF, see [11] or Appendix B. At this critical point the denominator of the RJ distribution vanishes exactly and the singularity of the RJ distribution is regularized by the macroscopic population of the fundamental mode. This leads to the following expression of the condensate fraction as a function of the energy:

$$ \frac{n_0^{eq}}{N} = 1 - \frac{E - E_0}{E_{crit}^* - E_0}. \qquad (46) $$

Then $n_0^{eq}$ vanishes at $E_{crit}^*$, and $n_0^{eq}/N \to 1$ as $E$ reaches the minimum value of the energy $E_0 = N\beta_0$, where $\beta_0$ is the fundamental mode eigenvalue, see Appendix B.

Because of finite size effects, the experiment reported in Ref.[12] does not occur in the strict thermodynamic limit. The theory of RJ condensation accounting for finite size effects gives the critical energy $E_{crit} = E_0(1 + (M - 1)/\varrho)$, where $\varrho$ is a constant number that tends to $V_0/\beta_0$ in the thermodynamic limit, see Appendix B. Considering the experimental parameters, we obtain $E_{crit}^*/E_{crit} \simeq 0.95$ (blue cross in Fig. 15b), so that the experiment is relatively 'close' to the thermodynamic limit.

In order to compare the theory of RJ condensation with the experiments, the condensate amplitude $n_0^{eq}$ and the chemical potential $\mu$ need to be extracted from the experimental data. In Ref.[12], this was done by implementing a fitting procedure. Figure 15a reports the chemical potential $\mu$ vs $E$: by decreasing $E$, $\mu$ increases and condensation occurs when $\mu \to \beta_0^-$ for $E = E_{crit}$. Below the transition ($E \le E_{crit}$) the fundamental mode gets macroscopically populated, $n_0^{eq} \gg n_{p \neq 0}^{eq}$, see Fig. 15b. The triangles report the experimental data retrieved
<figure><img src="image_16.png" /><figcaption>Figure 15: Rayleigh-Jeans condensation: (a) Chemical potential vs energy: For $E \lesssim E_{crit}$, $\mu \rightarrow \beta_0^-$, which leads to the macroscopic population of the fundamental mode, see $n_0^{eq}/N$ vs $E/E_{crit}$ (b). The blue (yellow) triangles report the experimental results from the NF (FF) intensity distributions (averaged over the realizations). The red lines report the RJ theory for the MMF used in the experiment. By decreasing the energy $E$ below the critical value $E_{crit}$, the chemical potential $\mu \rightarrow \beta_0^-$ (inset shows a zoom), which leads to the macroscopic population of the fundamental mode. The dashed black line in (b) refers to the thermodynamic limit (Eq.(46)) and the blue cross denotes $E_{crit}^*/E_{crit} \simeq 0.95$: The experiment is relatively 'close' to the thermodynamic limit. Source: From Ref.[12].</figcaption></figure>

from the fitting procedure from the near-field and far-field intensity distributions, see [12] for details. The red line reports the RJ theory accounting for finite size effects (i.e., beyond the thermodynamic limit) for the MMF used in the experiment, see Appendix B. The experimental results in Fig. 15 are in quantitative agreement with the RJ theory, where the parameters $\beta_0$ and $M$ are fixed by the fiber used in the experiment. Furthermore, the experimental results are close to the thermodynamic limit given by Eq.(46), see Fig. 15b. Notice that, as usual, finite size effects make the transition to condensation 'smoother', as evidenced by the red line in Fig. 15b.

### 3.3.2 Light thermalization to negative temperature equilibrium states

The experimental procedure implemented in Ref.[12] to study light thermalization with different energies has been extended to observe the thermalization to negative temperature equilibrium states. Indeed, equilibrium states with an energy larger than a threshold value $E_*(E_* = N \langle \beta_p \rangle)$, where $\langle \beta_p \rangle$ is the mean of the propagation constants, see Appendix B), are characterized by a negative temperature, i.e., an equilibrium state exhibiting an inverted modal population, where high-order modes are more populated than low-order modes. Originally conceived by Onsager [163], and Ramsey [164], negative temperatures created their own share of confusion for many decades, until being broadly accepted [165, 166], in line with recent experimental observations with cold atoms [167] and 2D quantum superfluids [168, 169]. More recently, negative temperatures have been predicted for classical multimode optical wave systems in the weakly nonlinear regime [13, 14], and subsequently observed in these systems with light waves in Refs.[38, 39]. Here, we summarize the experiments relevant to MMFs which have been reported in Ref.[38]. In these experiments, a GRIN MMF of 12 m was used, which guides nine groups of degenerate modes (with a
fiber radius of ≈ 15 µm). As described above in the experiments of RJ conden-
sation, an average over ensembles of speckle beams was considered: 2 × 300
realizations of the near-field and far-field intensity distributions were recorded
for the same power. The ensemble of realizations was then partitioned within
small energy intervals, and then averaged in each interval to get the modal
distributions reported in Fig. 16. The resolved modal distribution was obtained
from the well-known Gerchberg-Saxton algorithm [170, 171], which allows to
retrieve the transverse phase profile of the speckle field from the near-field and
far-field intensity distributions measured in the experiments, see Refs.[12, 141]
for details. The procedure was applied at the fiber output (12 m), as well as
the fiber input (after 20 cm of propagation). The results are reported in Fig. 16,
which evidence the process of light thermalization to negative temperature
equilibria that are characterized by an inverted modal population.

A distinguishing property of negative temperature equilibrium states is
that the radial near-field intensity distribution exhibits an oscillating behavior.
Indeed, the intensity distribution of usual positive-temperature equilibria is,
in average, a monotonic decreasing function with the radial distance. This is
consistent with the intuitive idea that low-order modes localized near the fiber
center are the most populated ones. In contrast, the inverted modal population
of negative temperature equilibria is characterized by an oscillating behavior
of the radial intensity, as illustrated in Fig. 17. Note that the number of radial
oscillations in Fig. 17 is given by the most oscillating mode of the fiber, that is
the mode LP₀₄ that exhibits five oscillations.

3.4 *Beyond thermodynamics: the modal phase-locking*

As discussed in the previous sections, the whole statistical mechanics ap-
proach to nonlinear multimode fiber systems relies on the pillar of averaging
over a meaningful statistical ensemble. As such, experiments of classical wave
thermalization and condensation involve results which are averaged over sev-
eral realizations. These have shed light on the description of the phenomenon
of BSC within the statistical mechanics framework, as discussed in Sec. 3.3.

At the same time, most of the experimental demonstrations of BSC did not
involve an explicit averaging process. Such an approach, which was used in
the experiments reported in Sec. 3.1 and Sec. 3.2, is justifiable by the fact
that averaging over several realizations is not practical for applications based
on BSC. Remarkably, however, the results of experiments carried out without
averaging over several realizations still lead to an excellent agreement with
some theoretical predictions, e.g., the match of the experimentally retrieved
mode occupancy with the RJ law. In this regard, a notable exception is the
entropy growth discussed in Sec. 3.1. Indeed, according to the statistical
mechanics theory, and specifically to the H-theorem, a monotonic evolution of
the entropy with the propagation (or the power) is expected when considering
an average over an ensemble of realizations of beam-cleaning experiments.
Whereas, the experimental results presented in Sec. 3.1 show that the trend of
the configuration entropy $\tilde{S}$ exhibits oscillations, i.e., it is not monotonic, at least
<figure><img src="image_17.png" /><figcaption>Figure 16: Thermalization to negative temperature equilibria. Experimental modal distributions averaged over realizations. Initial modal distributions at the fiber input (blue), and modal distributions at the fiber output (red). Corresponding RJ equilibrium distributions $n_p^{RJ}$ (green). The six panels correspond to different values of the energy (Hamiltonian) $E$, or equivalently different temperatures $T$ ($T > 0$ in (a), $T < 0$ in (b-f)). The modal distribution peaked on the lowest mode for $T > 0$ (a), gets inverted for $T < 0$ (b-f). An average over $\simeq 35$ realizations of speckle beams is considered for each panel. The fiber modes are sorted in the Hermite-Gauss basis from the fundamental one ($\beta_0$) to the highest mode group (nine-fold degenerate with $\beta_{max} = 9\beta_0$). Degenerate modes are equally populated at equilibrium, leading to a staircase distribution at equilibrium. Source: From Ref. [38].</figcaption></figure>

as far as out-of-equilibrium states are concerned (see Fig. 12). This seems to rule
out the hypothesis that the averaging operation is somehow intrinsically due to
the pulse-to-pulse variance of the mode-locked lasers used, for instance, in Ref.
[107] and Ref. [108]. Still, one may object that the oscillations of $\tilde{S}$ are not very
pronounced. In addition, the accuracy of holographic MD tools is questionable
when dealing with HOMs, which, at thermal equilibrium, provide the largest
contribution to the entropy.

In this regard, an important experimental proof that BSC, in the absence of
average over a statistical ensemble, cannot be fully identified as a wave ther-
malization phenomenon was presented in Ref. [111]. In that work, Mangini
et al. showed that BSC involves a mechanism of power-induced phase-locking
among the lowest-order modes, that carry most of the beam power. This ob-
<figure><img src="image_18.png" /><figcaption>Figure 17: Oscillating radial intensity distribution at negative temperature. Experimental averaged intensity distribution as a function of the radial (angle-averaged) distance |r| (red). Note the good agreement with the theoretical RJ intensity distribution (dashed green). The oscillating behavior of the intensity distribution is a signature of the negative temperature equilibrium. Inset: corresponding 2D intensity averaged over the realizations (the radius of the circle is the fiber radius). Source: From Ref.[38].</figcaption></figure>

servation collides with the statistical mechanics approach that only considers incoherent waves, i.e., it does not take into account the relative phase among modes.

Still, consistently with previous observations carried out with the same experimental setup [108], the mode occupancy distribution is in agreement with the statistical theory. This is shown in Fig. 18a,b, where the beam intensity profile and the mode power fraction ($A_j^2/P$, $A_j$ being the mode amplitude) at different input power levels are shown. As it can be seen, at low powers the output beam is speckled, meaning that thermal equilibrium was not yet reached. Whereas at high power levels the BSC effect takes place and the experimental mode amplitude agrees with the RJ law.

The establishment of RJ thermal equilibrium is nothing more than a mechanism of mode amplitude locking. Indeed, the amplitudes of the modes reach steady values as the input power grows larger. The left panel of Fig. 18c shows the amplitude locking effect. The equilibrium values of the amplitude are in agreement with the statistical mechanics predictions, as illustrated in the right panel of Fig. 18d. A similar locking effect for the mode amplitude was observed (cfr. left panel of Fig. 18d). Interestingly, by recurring to statistical mechanics considerations, it is possible to predict the values of the mode phases at thermal equilibrium. In particular, the experimental steady value of the mode phases turn out to be in excellent agreement with the predictions made with the following minimization problem

$$ \min_{\{\varphi_j\}} \{\varepsilon\} = \min_{\{\varphi_j\}} \{ |\Psi(\theta, r)|^2 - P \cdot |\psi_1(\theta, r)|^2 \}, \quad (47) $$

where $\Psi$ is the output intensity profiles of the total field and $\psi_1$ is fundamental mode intensity. As discussed in the previous sections, in fact, the beam temperature $T$ has to be sufficiently low, in order to ensure that a large amount of
the optical power is stored in the fundamental mode at thermal equilibrium. Hence, one can reasonably assume that at $T \approx 0^{+}$, the field $\Psi$ is close to the sole fundamental mode, which is the rationale behind the assumption of eq. (47).

However, it has to be mentioned that the predictive capability of the numerical minimization (47) turns out to be quite limited. In fact, the numerical procedure diverges whenever a relatively large number of modes is taken into account. Still, a unique solution can be found by limiting the problem to the six lowest-order modes (i.e., the first three mode groups), which, anyway, carry most of the power at thermal equilibrium in accordance to Fig. 18b. The remarkable effectiveness of eq. (47) in predicting the experimental outcomes can be visibly checked by comparing the left and right panel of Fig. 18d. In particular, it is remarkable to ascertain that modes belonging to the same group, which have almost the same amplitude, are nevertheless associated with different phases. As a final note, it is important to mention that an excellent agreement between experiments and numerical predictions was found for temperatures lower than $0.2 \text{ mm}^{-1}$, as illustrated in Fig. 18e.

In conclusion, BSC in GRIN MMFs results from a complex mode-locking phenomenon, where both mode amplitudes and spatial nonlinear phases remain simultaneously locked. The occurrence of a phase locking mechanism, in particular, is associated with the well-known spatial coherence preservation during BSC [85, 172]. This property is critically important for various applications, such as mode-locked fiber lasers and beam delivery in high-resolution nonlinear imaging, to name a few.

## 4 Conclusions and perspectives

In this review we presented a comprehensive analysis of optical beam propagation in nonlinear multimode fibers, through the lens of a statistical mechanics framework, with particular focus on the wave turbulence approach. This perspective has allowed us to elucidate the process of thermalization in GRIN MMFs, by capturing the intricate evolution of multimode laser beams through intermediate non-equilibrium states towards their final state of RJ equilibrium.

From the theoretical point of view, we reviewed recent works on the impact of random disorder, by extending wave turbulence theory to account for longitudinal refractive index fluctuations inherent to MMFs. This approach bridges the gap between theoretical predictions and the complex realities of light propagation in disordered environments, providing deeper insights into the role of environmental perturbations in beam dynamics.

As far as experimental investigations are concerned, we reviewed a number of key studies that have validated the statistical mechanics predictions for non-linear beam propagation. Specifically, we discussed recent direct observation of RJ mode distributions at positive, negative, and near-zero (wave condensation) temperatures. We also addressed the process of entropy growth during thermalization of either a single beam or two beams (optical calorimetry). In particular, the latter consists of the establishment of an equilibrium as a consequence of the heat exchange between two beams, which undoubtedly proves
<figure><img src="image_19.png" /><figcaption>Figure 18: (a) Experimental beam intensity profile when increasing *P*. (b) Comparison of the corresponding mode power distribution with the RJ law. (c,d) Comparison between experimental data (left) and predictions of the mode amplitude and phase, respectively. (e) Experimental (dots) and numerical (solid lines) values of the phase vs. beam temperature.</figcaption></figure>
that the thermodynamic parameters associated with the RJ mode distribution
have a true physical meaning.

Interestingly, while thermodynamics and wave turbulence models effec-
tively predict that mode amplitudes conform to the RJ law during BSC, statis-
tical mechanics arguments inherently fall short of capturing certain aspects of
self-cleaning, namely, the preservation of spatial coherence and the enhance-
ment of beam brightness. These, in fact, turn out to be associated with a modal
phase-locking mechanism, which can be remarkably well described by a simple
minimization problem (47). In this regard, we emphasize that the entropy does
not depend on the modal phases. Therefore, one may read Eq. (47) as being
analogous to a principle of minimum "energy" of a system with fixed entropy.
Moreover, BSC can be seen as the superposition of two thermodynamic effects
within two disjoint thermodynamic frameworks. The first involves the modal
amplitudes, and it relies on the maximization of entropy at fixed energy (the
Hamiltonian *H*). The second involves the modal phases, and it relies on the
minimization of an "energy", which is defined in (47) at fixed entropy. On the
other hand, we note that certain aspects of the WT theory have been recently
extended to account for the emergence of strongly correlated fluctuations in
turbulence flows [173, 166]. This could be an appropriate approach to develop
a detailed explanation of the effect of modal phase-locking observed experi-
mentally.

In terms of perspectives, many challenges await us. For instance, we foresee
that the next natural step will be the investigation of the statistical mechanics ap-
proach in the strongly nonlinear regime. Indeed, virtually all experiments that
demonstrated the process of RJ thermalization have been carried out in the weak
nonlinear regime, so far. Even a very recent experimental demonstration of the
Joule-Thomson expansion, which intrinsically deals with strong nonlinearities,
eventually relies on beam propagation in the weakly nonlinear regime [174].
The reason behind the researchers' fondness for the weakly nonlinear regime is
easy to guess: the presence of strong nonlinearity complicates the description
of the beam evolution toward thermal equilibrium, as well as the definition of
equilibrium state itself. Let us consider the case of wave condensation, which
is relatively simple to describe in cases where the nonlinearities are weak: in
the strong nonlinear regime, the condensed component forms a distinct phase,
characterized by different properties when compared to the thermal ( uncon-
densed) component. For instance, the presence of a strong condensate modifies
the linear dispersion relation of the thermal component, so that the condensate
amplitude n₀^(eq) is not simply given by the RJ distribution nₚ^(eq) taken at p = 0. One
needs to resort to a Bogoliubov-like transformation [175] in order to properly
describe the condensate fraction at equilibrium, see Ref.[57]. Notice that this
approach has been also extended to light propagation in GRIN MMFs, which
confirmed that the Bogoliubov correction is negligible for the weakly nonlinear
regime that it has been explored so far in optical thermalization experiments
[135]. It is also interesting to note that the modified Bogoliubov dispersion
relation has been recently measured experimentally in an optical system in a
bulk configuration [176]. Such a dispersion relation is at the origin of the pecu-
liar physical properties which have been experimentally investigated in fluids of light. Among these, we may cite superfluidity [177, 178, 179, 180], and the turbulence flow of superfluid light past an obstacle [44, 45], which are related to quantum vortex dynamics and the Berezinskii-Kosterlitz-Thouless transition [181].

As a future perspective, it would also be interesting to carry out optical ex-
periments aimed at studying the transition to condensation beyond the weakly
nonlinear regime explored so far in MMFs, in relation with the spontaneous
formation of strongly nonlinear coherent structures inherent to the focusing
regime (spatial solitons, rogue waves...), see e.g., [182].

Another aspect which has been largely disregarded so far and, we believe,
will be investigated in the near future, is the role that time plays in the non-
linear beam dynamics which leads to thermal equilibrium, e.g., in the case
of BSC. Indeed, since the results reported in Sec. 3.4 rule out the idea that
modal phases define the statistical ensemble that leads to the experimental
observation of the RJ distribution in standard BSC experiments, it is reason-
able to suppose that time averaging is the mechanism which is responsible for
the agreement between experimental observations and predictions of theories
based on statistical approaches. However, such hypothesis has not yet been
experimentally proved. As a matter of fact, only a few experimental works on
BSC have considered the time domain so far. Besides the results by Mangini et
al. discussed in Sec. 3, the experimental observations of Krupa et al. [183] and
Leventoux et al. [93], as well as Labaz and Sidorenko [184], indicate that the
temporal features of self-cleaned beam are non-trivial. Specifically, in [183, 93]
it was observed that spatial features of optical pulses that undergo BSC are
quite different when considering the center of the pulse or its tails; whereas
the Authors of [184] have shown that the spatial characteristics of self-cleaned
beams depend on the temporal profile of the input beam [184]. In addition,
it has to be noted that all experiments reported so far in the literature do not
take into account the intrinsic bandwidth of the measurement devices. Only
recently, Goery et al. presented at a conference a single-shot characterization
of the complex spatial amplitude of beams at the output of multimode fibers in
the nonlinear regime [185]. Unexpectedly, these results show that, when using
ultrashort pulses, the measurement bandwidth has a significant influence on
the output beam quality.

As far as temporal aspects are involved, it is worth emphasizing that the wave
turbulence theory discussed in Sec. 2 can be extended to study spatio-temporal
effects in MMFs [5, 186, 10, 87, 187, 188]. The development of a spatio-temporal
kinetic theory would also be important to study complex incoherent behaviors,
such as the generation of supercontinuum radiation in MMFs [5, 87, 189, 190,
191], in relation with the wave turbulence approach that has been developed
to study supercontinuum generation in singlemode fibers [50, 51, 53, 192]. The
discrete wave turbulence approach discussed in this work can also be extended
to study the impact of disorder on turbulence cascades, from both the theoretical
and the experimental point of view, in relation with the recent study in MMFs
[135].
Finally, a last research line that will undoubtedly be pursued is the im-
provement of the accuracy of MD techniques and the development of novel
MD devices. As a matter of fact, state-of-the-art MD tools have rather limited
accuracy when measuring the population of HOMs. This creates several issues.
A remarkable example in this regard is the divergence of the optical entropy
[88, 110]. Besides accuracy, future directions of MD involve the development of
fast devices and methods that, ideally, will allow to carry out the MD analysis
of single pulses as well as of time-resolved MD tools. In this regard, invaluable
help is likely to be provided by artificial intelligence algorithms, whose impact
on all fields of the scientific research has been exponentially growing over the
last few years.
Contents

<table><tbody><tr><td>1</td><td>Introduction and theoretical framework</td><td>2</td></tr><tr><td>1.1</td><td>The statistical mechanics theory</td><td>3</td></tr><tr><td>1.1.1</td><td>The wave turbulence theory</td><td>5</td></tr><tr><td>1.2</td><td>The beam self-cleaning effect</td><td>6</td></tr><tr><td>1.3</td><td>The goal and the content of this review</td><td>7</td></tr><tr><td>2</td><td>Wave turbulence in multimode fibers with disorder</td><td>9</td></tr><tr><td>2.1</td><td>Basic principles and conservation laws</td><td>9</td></tr><tr><td>2.2</td><td>Impact of random mode coupling</td><td>10</td></tr><tr><td>2.3</td><td>Kinetic equation with weak disorder</td><td>12</td></tr><tr><td>2.3.1</td><td>Model equation</td><td>12</td></tr><tr><td>2.3.2</td><td>Vanishing correlations among modes</td><td>15</td></tr><tr><td>2.3.3</td><td>Closure of the moments equations</td><td>15</td></tr><tr><td>2.3.4</td><td>Continuous wave turbulence</td><td>17</td></tr><tr><td>2.3.5</td><td>Discrete wave turbulence</td><td>19</td></tr><tr><td>2.3.6</td><td>Acceleration of thermalization mediated by disorder</td><td>21</td></tr><tr><td>2.3.7</td><td>Light thermalization resulting from turbulence cascades</td><td>22</td></tr><tr><td>2.3.8</td><td>Disorder-induced beam cleaning</td><td>25</td></tr><tr><td>2.3.9</td><td>Correlated and partially correlated disorder</td><td>26</td></tr><tr><td>2.3.10</td><td>Corrections of the dispersion relation</td><td>28</td></tr><tr><td>2.4</td><td>Kinetic equation with strong disorder</td><td>30</td></tr><tr><td>2.4.1</td><td>Model equation</td><td>30</td></tr><tr><td>2.4.2</td><td>Kinetic equation</td><td>31</td></tr><tr><td>2.4.3</td><td>Numerical simulations</td><td>32</td></tr><tr><td>2.5</td><td>Conclusions of the wave turbulence approach</td><td>36</td></tr><tr><td>3</td><td>Experimental evidences</td><td>38</td></tr><tr><td>3.1</td><td>The growth of optical entropy</td><td>40</td></tr><tr><td>3.1.1</td><td>Two-beam thermalization and optical calorimetry</td><td>42</td></tr><tr><td>3.2</td><td>High-temperature wave thermalization</td><td>45</td></tr><tr><td>3.3</td><td>Light thermalization and condensation by varying the beam energy</td><td>47</td></tr><tr><td>3.3.1</td><td>Rayleigh-Jeans condensation</td><td>47</td></tr><tr><td>3.3.2</td><td>Light thermalization to negative temperature equilibrium states</td><td>50</td></tr><tr><td>3.4</td><td>Beyond thermodynamics: the modal phase-locking</td><td>51</td></tr><tr><td>4</td><td>Conclusions and perspectives</td><td>54</td></tr><tr><td colspan="2">Appendix A Derivation of the discrete kinetic equation with weak disorder [Eq.(26)].</td><td>61</td></tr><tr><td colspan="2">Appendix A.1 Modal correlations</td><td>61</td></tr><tr><td colspan="2">Appendix A.2 Closure of the moments equations</td><td>63</td></tr><tr><td colspan="2">Appendix A.2.1 Computation of the term &lt;X<sub>p</sub><sup>(1)</sup>&gt;:</td><td>63</td></tr></tbody></table>
<table><tr><td>Appendix A.2.2</td><td>Computation of the term &lt;X<sub>p</sub><sup>(2)</sup>&gt;:</td><td>65</td></tr><tr><td>Appendix A.3</td><td>Impact of a correlated noise model of disorder on the kinetic equation</td><td>66</td></tr><tr><td>Appendix A.3.1</td><td>Computation of the moment &lt;(&#x2014;<sub>l</sub>&#x2014;<sub>m</sub>)(&#x2014;<sub>n</sub><sup>T</sup>&#x2014;<sub>p</sub>)&gt;</td><td>66</td></tr><tr><td>Appendix A.3.2</td><td>Computation of the moment &lt;(&#x2014;<sub>n</sub><sup>T</sup>&#x2014;<sub>m</sub>)(&#x2014;<sub>l</sub>&#x2014;<sub>p</sub>)&gt;</td><td>67</td></tr><tr><td>Appendix A.4</td><td>Impact of a partially correlated noise model of disorder on the kinetic equation</td><td>67</td></tr></table>

<table><tr><td>Appendix B</td><td>Wave condensation and condensate fraction plotted in Fig. 15.</td><td>69</td></tr><tr><td>Appendix B.1</td><td>Thermodynamic limit</td><td>69</td></tr><tr><td>Appendix B.2</td><td>Condensate fraction beyond the thermodynamic limit</td><td>70</td></tr></table>

<table><tr><td>Appendix C</td><td>Derivation of the kinetic equation with strong disorder [Eq.(37)].</td><td>71</td></tr><tr><td>Appendix C.1</td><td>Primary asymptotics</td><td>71</td></tr><tr><td>Appendix C.2</td><td>Secondary asymptotics</td><td>72</td></tr><tr><td>Appendix C.3</td><td>Degenerate modes</td><td>74</td></tr><tr><td>Appendix C.4</td><td>Model of disorder implemented in the simulations</td><td>76</td></tr></table>
# Appendix A Derivation of the discrete kinetic equation with weak disorder [Eq.(26)].

## Appendix A.1 Modal correlations

To study correlations among the modes, we derive an equation for the moments of the 2 × 2 matrix $\langle a_p^* a_q^T \rangle (z)$. We first compute the correlations within a specific mode ($p = q$), and next for distinct modes.

The second moments satisfy:

$$ \partial_z \langle a_p^* a_p^T \rangle = i \langle D_p^* a_p^* a_p^T \rangle - i \langle a_p^* a_p^T D_p^T \rangle - i \gamma \langle G_{pp}(a) \rangle. $$

where $G_{pq}(a(z)) = P_p(a)^* a_q^T(z) - a_p^* P_q(a)^T(z)$. According to the Furutsu-Novikov theorem:

$$ \langle \nu_{p,j} \sigma_j^* a_p^* a_p^T \rangle = \int_0^z \left\langle \frac{\delta(\sigma_j^* a_p^* a_p^T(z))}{\delta \nu_{p,j}(z')} \right\rangle \sigma_\beta^2 \mathcal{R}\left(\frac{z-z'}{l_\beta}\right) dz' $$

The variational derivative can be computed by following [193]. For $z > z'$ it is the solution of

$$ \partial_z \frac{\delta(a_p^* a_p^T(z))}{\delta \nu_{p,j}(z')} = i D_p^* \frac{\delta(a_p^* a_p^T(z))}{\delta \nu_{p,j}(z')} - i \frac{\delta(a_p^* a_p^T(z))}{\delta \nu_{p,j}(z')} D_p^T - i \gamma \frac{\delta(G_{pp}(a(z)))}{\delta \nu_{p,j}(z')} $$

starting from

$$ \frac{\delta(a_p^* a_p^T(z))}{\delta \nu_{p,j}(z')} \Big|_{z=z'} = i \sigma_j^* a_p^* a_p^T(z) - i a_p^* a_p^T(z) \sigma_j^T. $$

We need to know the form of the variational derivative for $z' < z$ and $|z - z'| = O(l_\beta)$. As $\sigma_\beta l_\beta \ll 1$ and $l_\beta \ll L_{nl}$, all terms in the right-hand side of the differential equation (A.1) satisfied by the variational derivative are negligible for $|z - z'| = O(l_\beta)$, so that the leading-order expression of the variational derivative for $z' < z$ and $|z - z'| = O(l_\beta)$ is

$$ \frac{\delta(a_p^* a_p^T(z))}{\delta \nu_{p,j}(z')} = i \sigma_j^* a_p^* a_p^T(z') - i a_p^* a_p^T(z') \sigma_j^T, \quad (A.1) $$

and therefore

$$ \langle \nu_{p,j} \sigma_j^* a_p^* a_p^T \rangle = i \int_0^z \langle \sigma_j^* \sigma_j^* a_p^* a_p^T(z') \rangle \sigma_\beta^2 \mathcal{R}\left(\frac{z-z'}{l_\beta}\right) dz' \\ - i \int_0^z \langle \sigma_j^* a_p^* a_p^T(z') \sigma_j^T \rangle \sigma_\beta^2 \mathcal{R}\left(\frac{z-z'}{l_\beta}\right) dz'. $$

For $j = 0$ this is zero and for $j = 1, 2, 3$ this can be approximated by (using $l_\beta \ll z$):

$$ \langle \nu_{p,j} \sigma_j^* a_p^* a_p^T \rangle = \frac{i \Delta \beta}{2} \left( \langle \sigma_j^* \sigma_j^* a_p^* a_p^T \rangle - \langle \sigma_j^* a_p^* a_p^T \rangle \langle \sigma_j^* \sigma_j^* \rangle \right), $$
and we find

$$ \partial_z \langle \mathbf{a}_p^* \mathbf{a}_p^T \rangle = -\Delta\beta \left( 3 \langle \mathbf{a}_p^* \mathbf{a}_p^T \rangle - \sum_{j=1}^3 \sigma_j^* \langle \mathbf{a}_p^* \mathbf{a}_p^T \rangle \sigma_j^T \right) - i\gamma \langle \mathbf{G}_{pp}(\mathbf{a}) \rangle. $$

The mean 2 × 2 matrix $\langle \mathbf{a}_p^* \mathbf{a}_p^T \rangle$ is Hermitian; therefore, it can be expanded as

$$ \langle \mathbf{a}_p^* \mathbf{a}_p^T \rangle = w_p(z) \sigma_0 + \sum_{j=1}^3 w_{p,j}(z) \sigma_j. $$

The real-valued functions $w_p$ and $(w_{p,j})_{j=1}^3$ satisfy

$$ \partial_z w_p = -\gamma \Pi_0 \{ \langle i \mathbf{P}_p(\mathbf{a})^* \mathbf{a}_p^T - i \mathbf{a}_p^* \mathbf{P}_p(\mathbf{a})^T \rangle \}, \quad (A.2) $$

$$ \partial_z w_{p,j} = -4\Delta\beta w_{p,j} - \gamma \Pi_j \{\langle i \mathbf{G}_{pp}(\mathbf{a}) \rangle\}, \quad (A.3) $$

where $\Pi_j W = \text{coefficient of the decomposition of the Hermitian matrix W on } \sigma_j$. In particular, $\Pi_0 W = \text{Tr}(W)/2$ so that

$$ \Pi_0 \{ \langle i \mathbf{G}_{pp}(\mathbf{a}) \rangle \} = -\text{Im} \langle \mathbf{P}_p(\mathbf{a})^\dagger \mathbf{a}_p \rangle. $$

The coefficients $w_{p,j}$ for $j = 1, 2, 3$ satisfy the damped equations

$$ \partial_z w_{p,j} = -4\Delta\beta w_{p,j} - \gamma \Pi_j \{\langle i \mathbf{G}_{pp}(\mathbf{a}) \rangle\}. $$

They are of the form

$$ w_{p,j}(z) = w_{p,j}(0) \exp(-4\Delta\beta z) - \gamma \int_0^z \exp(-4\Delta\beta(z-z')) \Pi_j \{\langle i \mathbf{G}_{pp}(\mathbf{a}) \rangle\} (z') dz'. $$

As we have $z, L_{nl} \gg 1/\Delta\beta$, the initial condition $w_{p,j}(0)$ is forgotten, and the second term in the right-hand side can be simplified and we get for $j = 1, 2, 3$

$$ w_{p,j} = -\frac{\gamma}{4\Delta\beta} \Pi_j \{\langle i \mathbf{G}_{pp}(\mathbf{a}(z)) \rangle\}. \quad (A.4) $$

Using the assumption $L_d = 1/\Delta\beta \ll L_{nl}$, we have

$$ w_{p,j} \approx 0 \quad (A.5) $$

to leading order.

Let us now study correlations among distinct modes ($p \neq q$). The variational derivatives satisfy for $z > z'$

$$ \partial_z \frac{\delta(\mathbf{a}_p^* \mathbf{a}_q^T(z))}{\delta \nu_{p,j}(z')} = i \mathbf{D}_p^* \frac{\delta(\mathbf{a}_p^* \mathbf{a}_q^T(z))}{\delta \nu_{p,j}(z')} - i \frac{\delta(\mathbf{a}_p^* \mathbf{a}_q^T(z))}{\delta \nu_{p,j}(z')} \mathbf{D}_q^T \\ + i(\beta_p - \beta_q) \frac{\delta(\mathbf{a}_p^* \mathbf{a}_q^T(z))}{\delta \nu_{p,j}(z')} - i\gamma \frac{\delta(\mathbf{G}_{pq}(\mathbf{a}(z)))}{\delta \nu_{p,j}(z')}, $$
with $G_{pq}(\mathbf{a}(z)) = P_p(\mathbf{a})^* \mathbf{a}_q^T(z) - \mathbf{a}_p^* P_q(\mathbf{a})^T(z)$ and starting from $\frac{\delta(\mathbf{a}_p^* \mathbf{a}_q^T(z))}{\delta\nu_{p,j}(z')} |_{z=z'} = \sigma_j^* \mathbf{a}_p^* \mathbf{a}_q^T(z')$. Assuming $\sigma_\beta l_\beta \ll 1$ and $l_\beta \ll L_{nl}$, the stochastic and nonlinear terms in the equation for the variational derivatives are negligible for $|z-z'| = O(l_\beta)$, so that to leading-order: $\delta(\mathbf{a}_p^* \mathbf{a}_q^T(z))/\delta\nu_{p,j}(z') = \sigma_j^* \mathbf{a}_p^* \mathbf{a}_q^T(z') \exp(i(\beta_p - \beta_q)(z-z'))$. Then using the involution property $\sigma_j^2 = \sigma_0$ we obtain $\langle D_p^* \mathbf{a}_p^* \mathbf{a}_q^T \rangle = 3i \int_0^z \langle \mathbf{a}_p^* \mathbf{a}_q^T(z') \rangle \sigma_\beta^2 \mathcal{R}\left(\frac{z-z'}{l_\beta}\right) \exp(i(\beta_p - \beta_q)(z-z')) dz'$, which can be approximated by

$$ \langle D_p^* \mathbf{a}_p^* \mathbf{a}_q^T \rangle = \frac{3i\Delta\beta}{2} \langle \mathbf{a}_p^* \mathbf{a}_q^T(z) \rangle, $$

by using $\sigma_\beta l_\beta \ll 1$ and $l_\beta \ll L_{nl}$ ($\mathbf{a}_p^* \mathbf{a}_q^T(z) \simeq \mathbf{a}_p^* \mathbf{a}_q^T(z') \exp(i(\beta_p - \beta_q)(z-z'))$ for $|z-z'| = O(l_\beta)$). Therefore we obtain

$$ \partial_z \langle \mathbf{a}_p^* \mathbf{a}_q^T \rangle = (-3\Delta\beta + i(\beta_p - \beta_q)) \langle \mathbf{a}_p^* \mathbf{a}_q^T \rangle - i\gamma \langle G_{pq} \rangle. $$

Writing the solution for $\langle \mathbf{a}_p^* \mathbf{a}_q^T \rangle (z)$ and using $z, L_{nl} \gg 1/\Delta\beta$, the initial condition $\langle \mathbf{a}_p^* \mathbf{a}_q^T \rangle (0)$ is forgotten and the solution can be written

$$ \langle \mathbf{a}_p^* \mathbf{a}_q^T \rangle (z) = \frac{i\gamma \langle G_{pq} \rangle}{3\Delta\beta - i(\beta_p - \beta_q)}. $$

Since $L_{nl} \gg 1/\Delta\beta$, we obtain $\langle \mathbf{a}_p^* \mathbf{a}_q^T \rangle (z) \simeq 0$.

## Appendix A.2 Closure of the moments equations

The diagonal modal components $w_p$ in (A.2) satisfy the undamped equation

$$ \partial_z w_p = \gamma \mathrm{Im} \langle \mathbf{P}_p(\mathbf{a})^\dagger \mathbf{a}_p \rangle. $$

With the expression (4) of the nonlinear term, we get that the mode occupancies $w_p$ satisfy the coupled equations (10-12). We now derive the equations governing the evolutions of the fourth-order moments $\langle X_p^{(1)} \rangle$ and $\langle X_p^{(2)} \rangle$ given by Eqs.(11-12).

### Appendix A.2.1 Computation of the term $\langle X_p^{(1)} \rangle$:

We first write the equation satisfied by the product of four vector fields:

$$ \begin{aligned} \partial_z (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) = & i(\beta_l + \beta_m - \beta_n - \beta_p)(\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) + i(\mathbf{a}_l^\dagger \mathbf{D}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \\ & +i(\mathbf{a}_l^\dagger \mathbf{D}_m^* \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) - i(\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{D}_n^T \mathbf{a}_p) \\ & -i(\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{D}_p \mathbf{a}_p) + i\gamma Y_{lmnp}^{(1)}, \end{aligned} \quad (A.6) $$
where

$$
\begin{equation}
\begin{aligned}
Y_{lmnp}^{(1)} = & - \sum_{l',m',n'} S_{l'm'n'l}^{*} \left[ \frac{1}{3} (\mathbf{a}_{l'}^{\dagger} \mathbf{a}_{m'}^{*}) (\mathbf{a}_{n'}^{T} \mathbf{a}_{m}^{*}) + \frac{2}{3} (\mathbf{a}_{n'}^{T} \mathbf{a}_{m'}^{*}) (\mathbf{a}_{l'}^{\dagger} \mathbf{a}_{m}^{*}) \right] (\mathbf{a}_{n}^{T} \mathbf{a}_{p}) \right. \\
& - \sum_{l',m',n'} S_{l'm'n'm}^{*} \left[ \frac{1}{3} (\mathbf{a}_{l'}^{\dagger} \mathbf{a}_{m'}^{*}) (\mathbf{a}_{l}^{\dagger} \mathbf{a}_{n'}) + \frac{2}{3} (\mathbf{a}_{n'}^{T} \mathbf{a}_{m'}^{*}) (\mathbf{a}_{l}^{\dagger} \mathbf{a}_{l'}) \right] (\mathbf{a}_{n}^{T} \mathbf{a}_{p}) \\
& + \sum_{l',m',n'} S_{l'm'n'n}^{*} \left[ \frac{1}{3} (\mathbf{a}_{l'}^{T} \mathbf{a}_{m'}) (\mathbf{a}_{n'}^{\dagger} \mathbf{a}_{p}) + \frac{2}{3} ((\mathbf{a}_{n'}^{\dagger} \mathbf{a}_{m'}) (\mathbf{a}_{l'}^{T} \mathbf{a}_{p})) \right] (\mathbf{a}_{l}^{\dagger} \mathbf{a}_{m}^{*}) \\
& \left. + \sum_{l',m',n'} S_{l'm'n'p}^{*} \left[ \frac{1}{3} (\mathbf{a}_{l'}^{T} \mathbf{a}_{m'}) (\mathbf{a}_{n}^{T} \mathbf{a}_{n'}^{*}) + \frac{2}{3} (\mathbf{a}_{n'}^{\dagger} \mathbf{a}_{m'}) (\mathbf{a}_{n}^{T} \mathbf{a}_{l'}) \right] (\mathbf{a}_{l}^{\dagger} \mathbf{a}_{m}^{*}). \right]
\end{aligned}
\tag{A.7}
\end{equation}
$$

We take the expectation and we apply the Gaussian summation rule to the
sixth-order moments in the expression of $\langle Y_{lmnp}^{(1)} \rangle$:

$$
\begin{align}
\langle Y_{lmnp}^{(1)} \rangle &= \frac{16}{3} S_{lmnp} (n_l n_m n_p + n_l n_m n_n - n_n n_p n_m - n_n n_p n_l) \nonumber \\
&\quad + \frac{16}{3} \delta_{mp}^K s_{ln}(\mathbf{w}) n_p (n_l - n_n) + \frac{16}{3} \delta_{mn}^K s_{lp}(\mathbf{w}) n_n (n_l - n_p) \nonumber \\
&\quad + \frac{16}{3} \delta_{lp}^K s_{mn}(\mathbf{w}) n_p (n_m - n_n) + \frac{16}{3} \delta_{ln}^K s_{mp}(\mathbf{w}) n_n (n_m - n_p), \tag{A.8} \\
s_{ln}(\mathbf{w}) &= \sum_{n'} S_{ln'n'n} n_{n'}. \tag{A.9}
\end{align}
$$

We now extend the procedure of Sec. Appendix A to the computation
of fourth-order modes considered here. Making use of the Furutsu-Novikov
theorem, and considering the general case $l \neq m \neq n \neq p$ in the regime
$\sigma_{\beta}l_{\beta} \ll 1$ and $l_{\beta} \ll L_{nl}$, we obtain

$$
\partial_z \langle (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \rangle = \left( -8\Delta\beta + i(\beta_l + \beta_m - \beta_n - \beta_p) \right) \langle (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \rangle \\
\qquad + i\gamma \langle Y_{lmnp}^{(1)} \rangle, \tag{A.10}
$$

whose solution has the form

$$
\begin{align*}
& \langle (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \rangle (z) \\
&= \langle (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \rangle (0) \exp \left( \left( -8\Delta\beta + i(\beta_l + \beta_m - \beta_n - \beta_p) \right) z \right) \\
&\quad + i\gamma \int_0^z \langle Y_{lmnp}^{(1)} \rangle (z') \exp \left( \left( -8\Delta\beta + i(\beta_l + \beta_m - \beta_n - \beta_p) \right) (z - z') \right) dz' \tag{A.11}
\end{align*}
$$

These equations correspond to those reported in Eq.(13) and Eq.(14) for the
fourth-order moment $J_{lmnp}^{(1)}$.

In the other cases when (at least) two indices are equal, i.e., the fourth-
order moments involve degenerate modes, the calculation shows that there is
still a damping in the moment equation, although the damping factor can be different from $8\Delta\beta$. However we will neglect this change because: (i) these terms are negligible in the triple sum and when the number of modes is large ($M \gg 1$), and (ii) this change only affects the multiplicative coefficient $8\Delta\beta$ in (A.10). Finally we give the expression of $\langle X_p^{(1)} \rangle$ in the discrete wave turbulence regime where $\beta_0 \gg \Delta\beta$. Collecting all terms we obtain

$$ \langle X_p^{(1)} \rangle = \frac{2\gamma}{3\Delta\beta} \sum_{l,m,n} |S_{lmnp}|^2 \delta^K(\beta_l + \beta_m - \beta_n - \beta_p) \\ \times (n_l n_m n_p + n_l n_m n_n - n_n n_p n_m - n_n n_p n_l) \\ + \frac{4\gamma}{3\Delta\beta} \sum_l |s_{lp}(\mathbf{w})|^2 \delta^K(\beta_l - \beta_p)(n_l - n_p). \quad (A.12) $$

## Appendix A.2.2 Computation of the term $\langle X_p^{(2)} \rangle$:

The equation satisfied by the product of four vector fields reads

$$ \partial_z (\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) = i(\beta_l + \beta_m - \beta_n - \beta_p)(\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) + i(\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{D}_l^\dagger \mathbf{a}_p) \\ +i(\mathbf{a}_n^T \mathbf{D}_m^* \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) - i(\mathbf{a}_n^T \mathbf{D}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) \\ -i(\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{D}_p \mathbf{a}_p) + i\gamma Y_{lmnp}^{(2)}, \quad (A.13) $$

where

$$ Y_{lmnp}^{(2)} = - \sum_{l',m',n'} S_{l'm'n'l}^* \left[ \frac{1}{3} (\mathbf{a}_l^\dagger \mathbf{a}_{m'}^*)(\mathbf{a}_{n'}^T \mathbf{a}_p) + \frac{2}{3} (\mathbf{a}_{n'}^T \mathbf{a}_{m'}^*)(\mathbf{a}_{l'}^\dagger \mathbf{a}_p) \right] (\mathbf{a}_m^\dagger \mathbf{a}_n) \\ - \sum_{l',m',n'} S_{l'm'n'm}^* \left[ \frac{1}{3} (\mathbf{a}_{l'}^\dagger \mathbf{a}_{m'}^*)(\mathbf{a}_{n'}^T \mathbf{a}_n) + \frac{2}{3} (\mathbf{a}_{n'}^T \mathbf{a}_{m'}^*)(\mathbf{a}_{l'}^\dagger \mathbf{a}_n) \right] (\mathbf{a}_l^\dagger \mathbf{a}_p) \\ + \sum_{l',m',n'} S_{l'm'n'n}^* \left[ \frac{1}{3} (\mathbf{a}_{l'}^T \mathbf{a}_{m'})(\mathbf{a}_{n'}^\dagger \mathbf{a}_m^*) + \frac{2}{3} (\mathbf{a}_{n'}^\dagger \mathbf{a}_{m'})(\mathbf{a}_{l'}^T \mathbf{a}_m^*) \right] (\mathbf{a}_l^\dagger \mathbf{a}_p) \\ + \sum_{l',m',n'} S_{l'm'n'p}^* \left[ \frac{1}{3} (\mathbf{a}_{l'}^T \mathbf{a}_{m'})(\mathbf{a}_{n'}^\dagger \mathbf{a}_l^*) + \frac{2}{3} (\mathbf{a}_{n'}^\dagger \mathbf{a}_{m'})(\mathbf{a}_{l'}^T \mathbf{a}_l^*) \right] (\mathbf{a}_m^\dagger \mathbf{a}_n). \quad (A.14) $$

We take the expectation and we apply the Gaussian summation rule to the sixth-order moments in the expression of $\langle Y_{lmnp}^{(2)} \rangle$:

$$ \langle Y_{lmnp}^{(2)} \rangle = \frac{16}{3} S_{lmnp} (n_l n_m n_p + n_l n_m n_n - n_n n_p n_m - n_n n_p n_l) \\ + \frac{16}{3} \delta_{mp}^K s_{ln}(\mathbf{w}) n_p (n_l - n_n) + \frac{32}{3} \delta_{mn}^K s_{lp}(\mathbf{w}) n_n (n_l - n_p) \\ + \frac{32}{3} \delta_{lp}^K s_{mn}(\mathbf{w}) n_p (n_m - n_n) + \frac{16}{3} \delta_{ln}^K s_{mp}(\mathbf{w}) n_n (n_m - n_p) \quad (A.15) $$

If $l \neq m \neq n \neq p$, then Furutsu-Novikov formula gives

$$ \partial_z \langle (\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) \rangle = (-8\Delta\beta + i(\beta_l + \beta_m - \beta_n - \beta_p)) \langle (\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) \rangle + i\gamma \langle Y_{lmnp}^{(2)} \rangle. $$
This equation corresponds to that reported in Eq.(13) for the fourth-order moment $J_{lmnp}^{(2)}$.

Following the same procedure, we have also derived the following equation

$$ \partial_z \langle (\mathbf{a}_l^T \mathbf{a}_l^*)(\mathbf{a}_m^\dagger \mathbf{a}_m) \rangle = 0, \qquad (A.16) $$

for any $l, m$. This implies, in particular, that the variance of the intensity fluctuations of each modal component $p$ is preserved $\langle |\mathbf{a}_p|^4 \rangle = \text{const}$, i.e., the Gaussian statistics is preserved during propagation.

Note that Eq.(A.16) is undamped, but this does not affect our results. Indeed, $\langle (\mathbf{a}_l^T \mathbf{a}_l^*)(\mathbf{a}_m^\dagger \mathbf{a}_m) \rangle$ is real-valued so that it does not contribute when it is substituted into Eq.(12) because $S_{lmml}$ is real-valued as well.

By neglecting the small corrections that appear in the cases when (at least) two indices are equal (i.e., fourth-order modes involving degenerate modes), we obtain in the discrete wave turbulence regime ($\beta_0 \gg \Delta\beta$):

$$ \langle X_p^{(2)} \rangle = \frac{2\gamma}{3\Delta\beta} \sum_{l,m,n} \delta^K(\beta_l + \beta_m - \beta_n - \beta_p) |S_{lmnp}|^2 \\ \times (n_l n_m n_p + n_l n_m n_n - n_n n_p n_m - n_n n_p n_l) \\ + \frac{2\gamma}{\Delta\beta} \sum_l \delta^K(\beta_l - \beta_p) |s_{lp}(\mathbf{w})|^2 (n_l - n_p). \qquad (A.17) $$

By replacing the expressions of $\langle X_p^{(j)} \rangle$ ($j = 1, 2$) given in (A.12) and (A.17) into the equation for the evolution of the modal components (10), we obtain the discrete kinetic Eq.(26).

## Appendix A.3 Impact of a correlated noise model of disorder on the kinetic equation

### Appendix A.3.1 Computation of the moment $\langle (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \rangle$

In the model of correlated disorder, all modes experience the same noise, $\mathbf{D}_n \equiv \mathbf{D} = \sum_{j=0}^3 \nu_j \sigma_j$. The equation (A.6) for the evolution of the product of four fields now reads

$$ \partial_z (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) = i(\beta_l + \beta_m - \beta_n - \beta_p)(\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) + i(\mathbf{a}_l^\dagger \mathbf{D}^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \\ +i(\mathbf{a}_l^\dagger \mathbf{D}^* \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) - i(\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{D}^T \mathbf{a}_p) \\ -i(\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{D} \mathbf{a}_p) + i\gamma Y_{lmnp}^{(1)}. \qquad (A.18) $$

We note that $\mathbf{a}_l^\dagger \mathbf{D}^\dagger \mathbf{a}_m^* + \mathbf{a}_l^\dagger \mathbf{D}^* \mathbf{a}_m^* = 2 \sum_{j \in \{1, 3\}} \nu_j \mathbf{a}_l^\dagger \sigma_j \mathbf{a}_m^*$. We follow the procedure outlined above in sections Appendix A.1-Appendix A.2. Using the Furutsu-Novikov theorem and assuming $\sigma_\beta l_\beta \ll 1$ and $l_\beta \ll L_{nl}$ we obtain

$$ \frac{\delta \mathbf{a}_l^\dagger \sigma_j \mathbf{a}_m^* \mathbf{a}_n^T \mathbf{a}_p(z)}{\delta \nu_j(z')} = 2i \left[ \mathbf{a}_l^\dagger \mathbf{a}_m^* \mathbf{a}_n^T \mathbf{a}_p(z') - \mathbf{a}_l^\dagger \sigma_j \mathbf{a}_m^* \mathbf{a}_n^T \sigma_j \mathbf{a}_p(z') \right] \\ \times \exp \left[ i(\beta_l + \beta_m - \beta_n - \beta_p)(z - z') \right] $$
and therefore

$$ \partial_z \langle (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \rangle = (-4\Delta\beta + i(\beta_l + \beta_m - \beta_n - \beta_p)) \langle (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \rangle \\ + i\gamma \langle Y_{lmnp}^{(1)} \rangle + 2\Delta\beta \sum_{j \in \{1,3\}} \langle (\mathbf{a}_l^\dagger \sigma_j \mathbf{a}_m^*)(\mathbf{a}_n^T \sigma_j \mathbf{a}_p) \rangle. \quad (A.19) $$

The last term can also be written as:

$$ \sum_{j \in \{1,3\}} \langle (\mathbf{a}_l^\dagger \sigma_j \mathbf{a}_m^*)(\mathbf{a}_n^T \sigma_j \mathbf{a}_p) \rangle = \langle (\mathbf{a}_l^\dagger \mathbf{a}_n)(\mathbf{a}_m^\dagger \mathbf{a}_p) \rangle + \langle (\mathbf{a}_l^\dagger \sigma_2 \mathbf{a}_n)(\mathbf{a}_m^\dagger \sigma_2 \mathbf{a}_p) \rangle. $$

The analysis reveals that the terms involving the dissipation that are proportional to $\Delta\beta$ in (A.19) essentially vanish. Indeed, using the factorizability property of statistical Gaussian fields to split the fourth-order moments into products of second-order moments, then we get $\langle (\mathbf{a}_l^\dagger \mathbf{a}_m^*)(\mathbf{a}_n^T \mathbf{a}_p) \rangle = \frac{1}{2} w_l w_m (\delta_{ln}^K \delta_{mp}^K + \delta_{lp}^K \delta_{mn}^K)$ and $\sum_{j \in \{1,3\}} \langle (\mathbf{a}_l^\dagger \sigma_j \mathbf{a}_m^*)(\mathbf{a}_n^T \sigma_j \mathbf{a}_p) \rangle = w_l w_m (\delta_{ln}^K \delta_{mp}^K + \delta_{lp}^K \delta_{mn}^K)$. This shows that, to leading order, the terms involving the dissipation indeed compensate with each other.

## Appendix A.3.2 Computation of the moment $\langle (\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) \rangle$

The equation (A.13) for the evolution of the product of four fields now reads

$$ \partial_z (\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) = i(\beta_l + \beta_m - \beta_n - \beta_p)(\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) + i(\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{D}^\dagger \mathbf{a}_p) \\ + i(\mathbf{a}_n^T \mathbf{D}^* \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) - i(\mathbf{a}_n^T \mathbf{D}^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) \\ - i(\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{D} \mathbf{a}_p) + i\gamma Y_{lmnp}^{(2)}. \quad (A.20) $$

Since the random matrix is Hermitian ($\mathbf{D}^\dagger = \mathbf{D}$), the expression can be reduced

$$ \partial_z (\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) = i(\beta_l + \beta_m - \beta_n - \beta_p)(\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) + i\gamma Y_{lmnp}^{(2)}. $$

The fourth-order moment $\langle (\mathbf{a}_n^T \mathbf{a}_m^*)(\mathbf{a}_l^\dagger \mathbf{a}_p) \rangle$ then evolves as in the absence of any disorder ($\Delta\beta = 0$).

## Appendix A.4 Impact of a partially correlated noise model of disorder on the kinetic equation

In the partially correlated model of disorder, all degenerate modes experience the same realization of disorder. Groups of degenerate modes with the same reduced eigenvalue $\beta$ are indexed by $p$, the modes within the $p$th group are indexed by $(p, j)$, and the linear polarization components of the $(p, j)$-th mode are $\mathbf{a}_{pj}$. The structural disorder induces a linear random coupling between the modes of different groups, as described by the $2 \times 2$ matrix $\mathbf{D}_p(z) = \sum_{l=0}^3 \nu_{p,l}(z) \sigma_l$ where $p$ labels the mode group number, i.e., degenerate modes that belong to the same group experience the same noise through the
random process $\nu_{p,l}(z)$ (this notation should not be confused with the decorrelated model of disorder where $p$ labels individual modes). The evolutions of the modal components are governed by

$$i\partial_z\mathbf{a}_{pj} = \beta_p\mathbf{a}_{pj} + \mathbf{D}_p(z)\mathbf{a}_{pj} - \gamma\mathbf{P}_{pj}(\mathbf{a}).$$

We look at the second-order moments $\langle \mathbf{a}_{pj}^* \mathbf{a}_{ql}^T \rangle$. We follow the procedure outlined above in section Appendix A.1 by using the Furutsu-Novikov theorem with $\sigma_\beta l_\beta \ll 1$ and $l_\beta \ll L_{nl}$.

Considering non-degenerate modes ($p \neq q$), we obtain

$$\langle \mathbf{a}_{pj}^* \mathbf{a}_{ql}^T \rangle (z) = \frac{i\gamma}{4\Delta\beta - i(\beta_p - \beta_q)} \langle \mathbf{P}_{pj}(\mathbf{a})^* \mathbf{a}_{ql}^T - \mathbf{a}_{pj}^* \mathbf{P}_{ql}(\mathbf{a})^T \rangle.$$

In the regime $L_d = 1/\Delta\beta \ll L_{nl}$ the correlation is vanishing small, $\langle \mathbf{a}_{pj}^* \mathbf{a}_{ql}^T \rangle (z) \approx 0$, as in the model of decorrelated disorder.

Let us now consider correlations among the orthogonal polarization components of a mode ($p = q$ and $j = l$), then we find as before

$$\langle \mathbf{a}_{pj}^* \mathbf{a}_{pj}^T \rangle (z) = w_{p,jl}(z)\sigma_0.$$

We consider correlations among distinct degenerate modes ($p = q$ and $j \neq l$).

We define the two Hermitian matrices $\mathbf{W}_{p,jl} = \frac{1}{2}(\langle \mathbf{a}_{pj}^* \mathbf{a}_{pl}^T \rangle + \langle \mathbf{a}_{pl}^* \mathbf{a}_{pj}^T \rangle)$ and $\tilde{\mathbf{W}}_{p,jl} = \frac{i}{2}(\langle \mathbf{a}_{pj}^* \mathbf{a}_{pl}^T \rangle - \langle \mathbf{a}_{pl}^* \mathbf{a}_{pj}^T \rangle)$ and then by carrying out the same calculations as in the case $p = q$ and $j = l$, we find

$$\langle \mathbf{a}_{pj}^* \mathbf{a}_{pl}^T \rangle (z) = w_{p,jl}(z)\sigma_0.$$

The coefficients $w_{p,jl}$ satisfy

$$\begin{aligned}
\partial_z w_{p,jl} &= \frac{\gamma}{2} \mathrm{Im} \langle \mathbf{P}_{pj}(\mathbf{a})^\dagger \mathbf{a}_{pl} + \mathbf{P}_{pl}(\mathbf{a})^\dagger \mathbf{a}_{pj} \rangle \mathbf{m} \\
&= \frac{1}{6} \gamma \langle X_{p,jl}^{(1)} + X_{p,lj}^{(1)} \rangle + \frac{1}{3} \gamma \langle X_{p,jl}^{(2)} + X_{p,lj}^{(2)} \rangle, \\
X_{p,jl}^{(1)} &= \mathrm{Im} \left\{ \sum_{p_1 j_1, p_2 j_2, p_3 j_3} S_{pl,p_1 j_1, p_2 j_2, p_3 j_3}^* (\mathbf{a}_{p_1 j_1}^\dagger \mathbf{a}_{p_2 j_2}^*) (\mathbf{a}_{p_3 j_3}^T \mathbf{a}_{pj}) \right\}, \\
X_{p,jl}^{(2)} &= \mathrm{Im} \left\{ \sum_{p_1 j_1, p_2 j_2, p_3 j_3} S_{pl,p_1 j_1, p_2 j_2, p_3 j_3}^* (\mathbf{a}_{p_3 j_3}^T \mathbf{a}_{p_2 j_2}^*) (\mathbf{a}_{p_1 j_1}^\dagger \mathbf{a}_{pj}) \right\}.
\end{aligned}$$

Let us look at the fourth-order moments $\langle (\mathbf{a}_{p_1 j_1}^\dagger \mathbf{a}_{p_2 j_2}^*) (\mathbf{a}_{p_3 j_3}^T \mathbf{a}_{p_4 j_4}) \rangle$ or $\langle (\mathbf{a}_{p_1 j_1}^T \mathbf{a}_{p_2 j_2}^*) (\mathbf{a}_{p_3 j_3}^\dagger \mathbf{a}_{p_4 j_4}) \rangle$. There are different types of such fourth-order moments that depend on the specific modes that they involve. Almost all of them satisfy an evolution equation with damping proportional to $\Delta\beta$, with different coefficients in front of $\Delta\beta$ that depend on the number of equal indices. These
terms are of the same form as those obtained by considering the decorrelated model of disorder, see section Appendix A.2. There are special cases when $p_j$ are equal by pairs (e.g. $p_1 = p_2$ and $p_3 = p_4$) where

$$ \partial_z \langle (\mathbf{a}_{p_1 j_1}^T \mathbf{a}_{p_1 j_2}^*) (\mathbf{a}_{p_3 j_3}^\dagger \mathbf{a}_{p_3 j_4}) \rangle = i \gamma Y_{p_3 j_3, p_1 j_2, p_1 j_1, p_3 j_4}^{(2)}, $$

which shows that there is no damping. From the expression of $Y^{(2)}$ (see Eq.(A.14)), and the evaluation of its expectation in terms of the $w_{p,jl}$ according to the Gaussian rule for sixth-order moments, such terms induce a reversible exchange of energy between the modes within each group.

## Appendix B Wave condensation and condensate fraction plotted in Fig. 15.

In this Appendix, we sketch the derivation of the condensation curve, i.e., the fraction of the power condensed in the fundamental mode vs the temperature, or the energy, in the classical RJ limit. This in relation with Section 3.3 where the experimental observation of RJ condensation has been reported.

### Appendix B.1 Thermodynamic limit

In order to render this appendix self-contained, we recall that we consider a 2D parabolic potential $V(\mathbf{r}) = q|\mathbf{r}|^2$ for $|\mathbf{r}| \le R$, that is truncated at $V_0 = qR^2$, where $R = 26\mu\text{m}$ is the fiber radius and $q$ a constant determined by the fiber characteristics, $q = k_0(n_{\text{co}}^2 - n_{\text{cl}}^2)/(2n_{\text{co}}R^2)$, $k_0 = 2\pi/\lambda_0$ being the laser wave-number, $n_{\text{co}} - n_{\text{cl}}$ the refractive index difference between the core and the cladding of the MMF. The eigenvalues are well approximated by the ones of the ideal harmonic potential $\beta_p = \beta_0(p_x + p_y + 1)$, which are solutions of the Schrödinger equation $\beta_p u_p(\mathbf{r}) = -\alpha\nabla^2 u_p(\mathbf{r}) + V(\mathbf{mr})u_p(\mathbf{r})$. The truncation of the potential $V(\mathbf{r}) \le V_0$ and the corresponding finite number of modes $M$ introduce an effective frequency cut-off in the far-field spectrum $k_c = \sqrt{2V_0/\beta_0}/r_0 \approx 1.15\mu\text{m}^{-1}$, where $r_0$ is the radius of the fundamental mode of the MMF [12].

We start from $N = \sum_p n_p^{eq}$ and $E = \sum_p \beta_p n_p^{eq}$ with the Rayleigh-Jeans (RJ) distribution $n_p^{eq} = T/(\beta_p - \mu)$, where $\beta_p = \beta_{p_x, p_y} = \beta_0(p_x + p_y + 1)$ are the eigenvalues of the truncated parabolic potential ($V(\mathbf{r}) \le V_0$), and the index $\{p\}$ labels the two integers $(p_x, p_y)$ that specify a mode. The sum over the modes reads $\sum_p = \sum_{0 \le p_x+p_y < g'}$, where $g = V_0/\beta_0$ is the number of groups of non-degenerate modes, with $M = g(g+1)/2$ the number of modes (see sec. III). Condensation arises when $\mu \to \beta_0$, which leads to a macroscopic population of the fundamental mode. In the following we term 'thermodynamic limit' the limit defined by $N \to \infty$, $\beta_0 \to 0$ with $N\beta_0^2 = \text{const}$ and $V_0 = \text{const}$. In this limit the discrete sums over the modes are replaced by continuous integrals, $N = \sum_p n_p^{eq} \to (T/\beta_0^2) \int_0^{V_0} dx \int_0^{V_0-x} dy (x+y+\bar{\mu})^{-1}$:

$$ N = \frac{T}{\beta_0^2} \left( V_0 + \bar{\mu} \log \left( \bar{\mu} / (V_0 + \bar{\mu}) \right) \right), \qquad (B.1) $$
with $\bar{\mu} = \beta_0 - \mu > 0$. Eq.(B.1) shows that $\mu \to \beta_0^-$ for a non-vanishing critical temperature

$$T_{\text{crit}}^{*} = N\beta_{0}^{2}/V_{0}. \qquad (B.2)$$

This means that wave condensation is reestablished in the thermodynamic limit in 2D owing to the truncated parabolic potential. The same conclusion is reached through the analysis of the energy. In the continuous limit we have

$$\frac{E}{N} = \mu + \frac{V_0^2/2}{V_0 + \bar{\mu} \log \left( \frac{\bar{\mu}}{V_0 + \bar{\mu}} \right)}. \qquad (B.3)$$

We see that $\mu \to \beta_0^-$ for a non-vanishing critical energy

$$E_{\text{crit}}^{*}/N = \beta_{0} + V_{0}/2 \approx V_{0}/2, \qquad (B.4)$$

and $\mu = \beta_0$ for $E \le E_{\text{crit}}^*$.

Note that, in contrast to homogeneous condensation (i.e. $V(\mathbf{r}) = 0$) where the density of states is constant, in the presence of the parabolic potential $V(\mathbf{r}) = q|\mathbf{r}|^2$ the density of states is a linear function of $\beta$ ($\rho(\beta) = \beta/\beta_0^2$), which prevents the infrared divergence of the integral of the power $N = \int_0^{V_0} \rho(\beta)n_\beta^{eq}d\beta$.

## Appendix B.2 Condensate fraction beyond the thermodynamic limit

The experiments are not performed in the thermodynamic limit, and finite size effects are taken into account by considering discrete sums over the modes (instead of continuous integrals): $N = \sum_p n_p^{eq}$, $E = \sum_p \beta_p n_p^{eq}$. As in the usual treatment of condensation, we split the contribution of the fundamental mode:

$$N = n_0^{eq} + T \sum_{p \neq 0} 1/(\beta_p - \beta_0), \qquad (B.5)$$

$$E = n_0^{eq} \beta_0 + T \sum_{p \neq 0} \beta_p / (\beta_p - \beta_0). \qquad (B.6)$$

By taking the ratio we eliminate the temperature $T$, which gives

$$\frac{n_0^{eq}}{N} = 1 - \frac{\varrho(E - N\beta_0)}{N\beta_0(M - 1)}, \qquad (B.7)$$

where $\varrho = \sum_{p \neq 0} 1/(p_x + p_y)$. Eq.(B.7) shows that condensation arises by decreasing the energy below the critical value

$$E_{\text{crit}} = N\beta_0(1 + (M - 1)/\varrho). \qquad (B.8)$$

Following the same procedure with the temperature instead of the energy, we obtain

$$n_0^{eq}/N = 1 - T/T_{\text{crit}}, \quad T_{\text{crit}} = N\beta_0/\varrho. \qquad (B.9)$$