# Rayleigh-Jeans thermalization vs beam cleaning in multimode optical fibers

K. Baudin, J. Garnier, A. Fusaro, K. Krupa, C. Michel, G. Millot, Antonio

Picozzi

► To cite this version:

K. Baudin, J. Garnier, A. Fusaro, K. Krupa, C. Michel, et al. Rayleigh-Jeans thermalization vs beam cleaning in multimode optical fibers. Optics Communications, 2023, 545, pp.129716. 10.1016/j.optcom.2023.129716. hal-04198185

HAL Id: hal-04198185

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

Submitted on 6 Sep 2023

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>
# Rayleigh-Jeans thermalization vs beam cleaning in multimode optical fibers

K. Baudin¹, J. Garnier², A. Fusaro³, C. Michel⁴, K. Krupa⁵, G. Millot¹,⁶, A. Picozzi¹

¹ Laboratoire Interdisciplinaire Carnot de Bourgogne, CNRS, Université de Bourgogne, Dijon, France

² CMAP, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, 91128 Palaiseau Cedex, France

³ CEA, DAM, DIF, F-91297 Arpajon Cedex, France

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

⁵ Institute of Physical Chemistry Polish Academy of Sciences, Warsaw, Poland and

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

Classical nonlinear waves exhibit, as a general rule, an irreversible process of thermalization toward the Rayleigh-Jeans equilibrium distribution. On the other hand, several recent experiments revealed a remarkable effect of spatial organization of an optical beam that propagates through a graded-index multimode optical fiber (MMF), a phenomenon termed beam self-cleaning. Our aim here is to evidence the qualitative impact of disorder (weak random mode coupling) on the process of Rayleigh-Jeans thermalization by considering two different experimental configurations. In a first experiment, we launch speckle beams in a relatively long MMF. Our results report a clear and definite experimental demonstration of Rayleigh-Jeans thermalization through light propagation in MMFs, over a broad range of kinetic energy (i.e., degree of spatial coherence) of the injected speckle beam. In particular, the property of energy equipartition among the modes is clearly observed in the condensed regime. The experimental results also evidence the double turbulence cascade process: while the power flows toward the fundamental mode (inverse cascade), the energy flows toward the higher-order modes (direct cascade). In a 2nd experiment, a coherent laser beam is launched into a relatively short MMF length. It reveals an effect of beam cleaning driven by an incipient process of Rayleigh-Jeans thermalization. As discussed through numerical simulations, the fast process of Rayleigh-Jeans thermalization observed in the 1st experiment can be attributed due to a random phase dynamics among the modes, which is favoured by the injection of a speckle beam and the increased impact of disorder in the long fiber system.

## I. INTRODUCTION

Understanding the mechanisms responsible for self-organization processes in conservative wave systems is a difficult problem that generated significant interest. Contrary to a dissipative system, a conservative Hamiltonian system cannot evolve towards a fully ordered state, because such an evolution would imply a loss of statistical information that would violate the formal reversibility of the evolution equations. However, in spite of its formal reversibility, a nonintegrable Hamiltonian wave system is expected to exhibit an irreversible evolution toward the equilibrium state, i.e., the most disordered state that realizes the maximum of entropy [1]. The *thermalization* to the Rayleigh-Jeans (RJ) equilibrium distribution of weakly nonlinear random waves can be described in detail by the well-developed wave turbulence theory [2, 3], that has been also successfully applied to optical waves [4–6].

RJ thermalization can be characterized by a counter-intuitive process of self-organization. To discuss this aspect with some concrete examples, let us consider the nonlinear Schrödinger equation (NLS). In the focusing regime, the self-organization process is known as “soliton turbulence”. It is characterized by the spontaneous formation of a large scale coherent solitary-wave that remains immersed in a sea of small-scale fluctuations [3, 4, 7–15]. On the other hand, in the defocusing regime of the NLS equation, the self-organization process refers to ‘wave condensation’ [3, 4, 16–27]. It originates in the

singular behaviour of the RJ equilibrium distribution, in analogy with the quantum Bose-Einstein condensation. More precisely, the RJ distribution exhibits a phase transition to a condensed state that is characterized by a macroscopic population of the fundamental mode of the system.

The observation of thermalization with optical waves in a (cavityless) free propagation is known to require very large propagation lengths, a feature that has been discussed recently in different circumstances [25, 28]. The situation is completely different when the optical beam propagates in a waveguide configuration. Indeed, the finite number of modes supported by the waveguide significantly increases the rate of thermalization to the RJ distribution [29]. In addition, the finite number of modes introduces a frequency cut-off that regularizes the ultraviolet catastrophe inherent to classical waves [29, 30]. It is in this framework that a remarkable phenomenon of spatial beam self-organization, called ‘beam self-cleaning’, has been recently discovered in graded-index multimode fibers (MMFs) [31–34]. Recent works indicate that this phenomenon of beam cleaning can be interpreted as a consequence of RJ thermalization and condensation, a feature that has been discussed both theoretically [35–39] and experimentally [40–47]. We recall in this respect that light propagation in an optical fiber is affected by a structural disorder due to refractive index fluctuations introduced by inherent imperfections and environmental perturbations [48–51, 53, 54], which significantly impact the process of RJ thermalization [35, 46, 52]. We also note that RJ thermalization to negative tempera-
ture equilibrium states in multimode systems has been theoretically predicted [38, 43], and recently observed in graded-index multimode optical fibers [55].

In this work we report an accurate experimental demonstration of the process of RJ thermalization through light propagation in a MMF. In a 1st experiment, we consider the injection of speckle beams in a relatively long MMF [40]. We compare the experimental distribution of the power among the fiber modes with the corresponding theoretical RJ distribution, which reveals a quantitative agreement without using adjustable parameters. Furthermore, such a remarkable agreement is obtained over a broad range of kinetic energies $E$ ('temperatures'), i.e., over a broad range of randomness of the speckle beams launched in the MMF. In this way, we provide the experimental observation of the equipartition of the kinetic energy among the fiber modes. The experimental results also evidence that, while the power (or 'wave-action', or particle number in a corpuscular picture) flows toward the fundamental mode (inverse turbulence cascade), the energy flows toward the high-order modes (direct turbulence cascade). Next, in a 2nd experiment, we consider the injection of a coherent laser beam in a short piece of MMF length, while the power is increased so as to keep constant the effective nonlinear interaction with respect to the 1st experiment. At variance with the 1st experiment, in the 2nd experimental configuration there are no random phases among the initially populated modes, while the impact of disorder on light propagation (polarization mixing and random mode coupling) is severely limited by the short fiber length used. In this 2nd experiment, while we can identify an effect of beam self-cleaning driven by the thermalization of the low-order modes, the high-order modes exhibit a deviation from the RJ distribution. As will be discussed through numerical simulations, the fast process of Rayleigh-Jeans thermalization observed in the 1st experiment can be attributed due to a random phase dynamics among the modes, which is favoured by the injection of a speckle beam and the increased impact of disorder in the long fiber system.

## II. FIRST EXPERIMENT: RAYLEIGH-JEANS THERMALIZATION

### A. Experimental setup and procedure

Figure 1 schematically reports the experimental setup. The source is a Nd:YAG laser delivering subnanosecond pulses (400ps) at $\lambda_0 = 1064$ nm. A key feature of the 1st experiment is that the laser beam is passed through a diffuser before the injection of the speckle beam into the MMF. More specifically, the laser beam is collimated and passed through a glass diffuser plate placed in the vicinity of the Fourier plane of a 4f-optical system. The near-field (NF) intensity distribution of the output beam was magnified and imaged on a first CCD camera owing

<figure><img src="image_2.png" /><figcaption>FIG. 1: **Experimental setup:** laser, optical isolator, half-wave plate and polarizer, diffuser (D), lenses for magnification and imaging ($f_j$), graded-index MMF, and cameras for near- and far-field detections (C).</figcaption></figure>

to a two lens telescope optical system. The CCD camera was placed on a rail orthogonal to the beam propagation in order to remove or put the camera back on the beam path. The far-field (FF) intensity distribution of such a magnified image was obtained by placing it in the object focal-plan of a lens and using a second CCD camera positioned in its image (Fourier) focal-plan, see Fig. 1.

An other important property of the 1st experiment considered in this section is that the MMF used is relatively long, i.e., $L = 12$ m long graded-index MMF. Its refractive index profile exhibits a parabolic shape in the fiber core with a maximum core index (at the center) of $n_{co}=1.470$ and $n_{cl} = 1.457$ for the cladding at the pump wavelength of 1064nm (numerical aperture NA=0.195, fiber radius $R = 26\mu$m). The corresponding parabolic potential reads $V(r) = q|r|^2$ for $|r| \le R$, where $R = 26\mu$m is the fiber radius and $q = k_0(n_{co}^2 - n_{cl}^2)/(2n_{co}R^2)$, $k_0 = 2\pi/\lambda_0$ the laser wave-number. The MMF guides $M \simeq 120$ modes (i.e., $g_{max} = 15$ groups of non-degenerate modes). The truncation of the potential introduces a frequency cut-off in the FF spectrum $k_c = (2\pi/\lambda_0)\sqrt{n_{co}^2 - n_{cl}^2}$. The eigenvalues are well approximated by the ones of the ideal harmonic potential $\beta_p = \beta_0(p_x+p_y+1)$, where the index $\{p\}$ labels the two numbers $(p_x, p_y)$ needed to specify a mode. We have computed analytically the propagation of the optical wave throughout the setup of our detection scheme. The optical field amplitude in the NF plane is an exact magnification of the wave amplitude at the output of the MMF. In addition, the optical field amplitude in the FF detection plane exactly corresponds to the Fourier transform of the amplitude at the fiber output. We note that the experimental setup for the NF and FF detection does not introduce detrimental spurious transverse phase profiles in the optical field, e.g., related to optical free propagation in air or phase shifts due to the presence of additional lenses. For more details about the experimental set-up, we refer the reader to Ref. [55] (Supplementary Material).

This setup allowed us to retrieve an accurate measurement of the power $N$ and the energy $E$ of the speckle beam. More precisely, the NF intensity distribution $I_{NF}(r) = |\psi|^2(r)$ provides a measurement of the power $N = \int I_{NF}(r)dr$ and of the potential energy $E_{pot} = \int V(r)|\psi|^2(r)dr$. The kinetic energy $E_{kin} = \alpha \int |\nabla\psi|^2(r)dr$ is retrieved from the FF intensity distribution $I_{FF}(k) = |\tilde{\psi}|^2(k)$, where $\alpha = 1/(2n_{co}k_0)$ and $\tilde{\psi}(k)$ is the Fourier transform of $\psi(r)$. This provides
<figure><img src="image_3.png" /><figcaption>FIG. 2: 1st Experiment: Rayleigh-Jeans thermalization. Modal distribution $\tilde{n}_g$ vs mode group g in normal scale (1st column), in log-scale (2nd column), and corresponding energy distribution $\tilde{E}_g$ vs g (3rd column). Initial modal distribution at the fiber input (blue), and modal distributions at the fiber output (red). Corresponding RJ distribution (green line): The quantitative agreement with the output distribution is obtained without using any adjustable parameter. Different lines in the figure correspond to different energies E, i.e., different condensate fractions (see the legend in the 1st column).</figcaption></figure>

the measurement of the (linear) energy (Hamiltonian) $E = E_{\text{pot}} + E_{\text{kin}}$. Projecting onto the basis of the fiber modes, the power and (kinetic) energy read

$$N = \sum_{p} n_{p}, \quad E = \sum_{p} \beta_{p} n_{p}, \qquad (1)$$

where $n_p$ is the amount of power in the mode $\{p\}$, and $\sum_p$ is carried over the set of M modes indexed by $\{p\}$ [29]. Note that E is in units of W.m$^{-1}$ (not Joule). However, we refer E to the energy because it is the (dominant) linear contribution to the Hamiltonian.

The accurate measurements of the NF and FF intensity distributions allowed us to retrieve the modal content of the optical beam in the experiment, namely the distribution of power among the modes $n_p^{\text{exp}}$. In general, to achieve the mode decomposition, several interferometric approaches based on use of a reference beam have been exploited to study light thermalization in MMFs [42, 44]. On the other hand, a number of methods without a reference

beam have been also proposed, see references in [56]. Here, we use a (non-interferometric) numerical computing mode decomposition procedure that is based on the well-established Gerchberg-Saxton algorithm [57, 58]. It allows us to retrieve the transverse phase profile of the field from the NF and the FF intensity distributions measured in the experiments [59, 60]. By projecting the retrieved spatial profile of the complex field over the fiber modes, we get the complete modal distribution of the experimental optical beam $n_p^{\text{exp}}$. We refer the reader to the Ref.[55] (Supplementary Material) for a detailed presentation of the procedure (theoretical and numerical computation of the error introduced by the Gerchberg-Saxton algorithm, impact of the sampling due to the cameras, validity of the phase reconstruction). In addition, we show experimentally and theoretically in the Appendix that the modal distribution recorded experimentally $n_p^{\text{exp}}$ at the fiber output converges toward the expected RJ distribution as the number of realizations of speckles beams is increased.
<figure><img src="image_4.png" /><figcaption>FIG. 3: 1st Experiment: Attraction to the RJ distribution. (a) Power fraction into the fundamental mode (condensate fraction) $n_0/N$ vs $E/E_{crit}$: at the fiber input (blue), and fiber output (red). Corresponding RJ theory (green line): The quantitative agreement with the output distribution is obtained without using any adjustable parameter. (b) Distance $D$ vs energy [defined from Eq.(4)] to the RJ equilibrium distribution at the fiber input (blue), and fiber output (red). The strong reduction of $D$ from the input-to-output shows the attraction to the RJ distribution. Note that $D$ is normalized, $0 \le D \le 1$.</figcaption></figure>

## B. RJ thermalization over a broad energy range

The experiment is realized in the weakly nonlinear regime, where linear propagation effects dominate over nonlinear effects $L_{lin} \sim \beta_0^{-1} \ll L_{nl} = 1/(\gamma N)$, where $\gamma$ is the nonlinear coefficient of the MMF. Typically, we have $L_{lin} \sim 0.2$mm, and $L_{nl} \sim 0.3$mm. This means that the nonlinearity plays a role over a propagation length much larger than that of the rapidly fluctuating speckle beam.

The process of RJ thermalization in the MMF is driven by the Kerr nonlinearity, as described theoretically by a wave turbulence kinetic approach [29, 35, 37, 46]. Accordingly, as the beam propagates through the MMF, it is expected to relax toward the thermodynamic equilibrium state described by the RJ distribution [29, 35, 37, 38, 40, 42–44].

$$n_p^{\text{RJ}} = T/(\beta_p - \mu), \quad (2)$$

where $T$ and $\mu$ are the temperature and chemical potential. Accordingly, we have at equilibrium $N = T \sum_p (\beta_p - \mu)^{-1}$ and $E = T \sum_p \beta_p / (\beta_p - \mu)$. Important to notice, $(T, \mu)$ are uniquely determined by $(N, E)$, i.e., there is a one to one correspondence between $(T, \mu)$ and $(N, E)$. In other terms, we are dealing with the microcanonic statistical description, so that $T$ does not refer to the 'room' temperature ($T$ is in units of W.m$^{-1}$). Also note that the RJ distribution can be generalized by taking into account for the conservation of linear momentum [61–64] and angular momentum [45, 65]. However, the conservation of angular momentum is relevant when peculiar conditions of injection of the optical beam into the MMF are considered, see Ref.[45]. Here, we consider speckle beams for which the angular momentum may vary from one realization to another and vanishes in average, so that its conservation is not relevant to our experiments.

At variance with previous experiments of spatial beam

cleaning and RJ thermalization [31, 33, 34, 40–43], here we study the process of thermalization for different amount of the (kinetic) energy $E$. Indeed, by passing the laser beam through a diffuser before injection in the MMF, we can vary the amount of randomness (i.e., degree of spatial coherence) of the speckle beam by keeping fixed the power $N$ – the higher the randomness of the fluctuations in the speckle beam, the higher the corresponding energy $E$. In other terms, owing to the diffuser, we study the process of RJ thermalization over a broad range of energies $E$. This is at variance with previous experiments [31, 33, 34, 40–43], as will be discussed in the framework of the 2nd experiment presented below in the paper. Note that throughout the paper, the terminology 'coherence' refers to the purely spatial coherence properties of the quasi-monochromatic optical beam that propagates through the fiber.

An other important advantage of using a diffuser is that it allows us to perform an average over an ensemble of realizations. Indeed, the RJ distribution is in essence a statistical distribution, so that the comparison with the experimental data requires an average over realizations. We have recorded 2×1000 realizations of the NF and FF intensity distributions for the same power $N$ ($N = 7$ kW) and different values of the energy $E$. For each individual realization, we have retrieved the modal distribution $|a_p^{\text{exp}}|^2$. We have partitioned such an ensemble of 1000 realizations of $\{|a_p^{\text{exp}}|^2\}$ within small intervals of energies $\Delta E = 0.125 E_0$, where $E_0 = N\beta_0$ is the minimum of the energy (only the fundamental mode is populated). Next, we have performed an average over the realizations of the modal distributions for each specific energy interval $\Delta E$, which gives the averaged modal distributions $n_p^{\text{exp}}$ for different energies $E$. We have followed this experimental procedure at the output of the MMF (at $L = 12$m), and then we have repeated the procedure at the fiber input (after 20cm of propagation).

### 1. Power distribution among the modes

We report in Fig. 2 (left column) the averaged modal distribution in the experiment $n_p^{\text{exp}}$ at the fiber input and the fiber output, for different energies $E$. For convenience we have reported the average power $\tilde{n}_g$ within each group of degenerate modes, where $g = 0, \dots, g_{max} - 1$ indexes the mode group. For the MMF used in the experiments, $g_{max} = 15$ for $M = g_{max}(g_{max} + 1)/2 = 120$ modes.

We compare the experimental data with the theoretical RJ distribution $n_p^{\text{RJ}}$. Note that there are no adjustable parameters between the experimental distribution $n_p^{\text{exp}}$ and the RJ distribution $n_p^{\text{RJ}}$: The parameters $(T, \mu)$ in $n_p^{\text{RJ}}$ in Eq.(2) are uniquely determined by the power and the energy $(N, E)$ measured in the experiments. We observe in Fig. 2 a very good agreement between the observed modal populations $n_p^{\text{exp}}$ and the RJ distribution $n_p^{\text{RJ}}$. Furthermore, such a good agreement is obtained
over a broad range of the energy $E$.

In order to properly appreciate the key role of the energy $E$, we report in Fig. 3(a) the condensate fraction (i.e., amount of power fraction in the fundamental fiber mode) $n_0^{\text{RJ}}/N$ vs $E/E_{\text{crit}}$, where $E_{\text{crit}} \simeq N\beta_0\sqrt{M/2}$ is the threshold value of the energy below which the fundamental mode gets macroscopically populated [40]. In a loose sense, the energy $E$ plays a role analogous to the temperature, i.e., the condensate fraction $n_0^{\text{RJ}}/N$ increases as the energy $E$ decreases. Note in particular that the effect of beam self-cleaning is suppressed nearby the transition to condensation for $E \simeq E_{\text{crit}}$, a feature discussed in detail in Ref. [40].

## 2. Energy distribution among the modes

The modal distribution of the power is mainly sensitive to the behaviour of the low-order modes. So as to properly characterize the thermalization process over all modes, it proves convenient to consider the distribution of the kinetic energy $E_p = (\beta_p - \beta_0)n_p$, or equivalently the average energy within groups of degenerate modes $\tilde{E}_g = \beta_0 g \tilde{n}_g$, where we recall that $g$ indexes the mode group. The results reported in Fig. 2 (right column) confirm a remarkable quantitative agreement between the experimental data and the RJ equilibrium distribution. We note in particular that for a small value of the energy $E \ll E_{\text{crit}}$, we can observe an energy equipartition among the modes, a property inherent to classical statistical mechanics. Indeed, well below the energy threshold for condensation, i.e. $E \ll E_{\text{crit}}$, the chemical potential reaches the fundamental mode eigenvalue $\mu \to \beta_0^-$, so that $n_p^{\text{RJ}} = T/(\beta_p - \beta_0)$ for $p \neq 0$ (or $\tilde{n}_g = T/(g\beta_0)$ for $g \neq 0$), and thus

$$E_p = T \quad \text{for } p \neq 0 \qquad (3)$$

(or $\tilde{E}_g = T$ for $g \neq 0$). This property of energy equipartition can be observed in all panels (k-n), $\tilde{E}_g^{\text{exp}} = T$, the only panel which shows a significant departure from energy equipartition is panel (o) because $E$ is then close to $E_{\text{crit}}$. Note that this is consistent with the fact that for $E \gtrsim E_{\text{crit}}$, the chemical potential deviates from the fundamental mode eigenvalue, and the RJ distribution (2) with $\mu \ll \beta_0$ no longer verifies the property of energy equipartition among the modes.

The experimental results also evidence the double turbulence cascade underlying RJ thermalization: While the power flows toward the fundamental mode (inverse cascade), the flow of kinetic energy toward the higher-order modes (direct cascade). This latter effect is shown in panels (l-o) of Fig. 2, for a relative small condensate fraction $n_0^{\text{RJ}}/N = 0.15$. This double cascade process was already commented in Ref. [43] through the analysis of the experimental data in direct $r$-space. The mode-resolved analysis reported in Fig. 2 is more appropriate and convincing.

<figure><img src="image_5.png" /><figcaption>FIG. 4: 2nd Experiment: Incipient RJ thermalization. (a) Distance $D$ vs energy $E/E_{\text{crit}}$ [defined from Eq. (4)] to the RJ equilibrium distribution at the fiber input (blue), and fiber output (red). (b)-(d) Modal distribution $\tilde{n}_g$ vs mode group $g$ (right column), and corresponding energy distribution $\tilde{E}_g$ vs $g$. Initial modal distribution at the fiber input (blue), and modal distributions at the fiber output (red), and corresponding RJ distribution (green line). Note that $D$ in (a) is bounded, $0 \le D \le 1$.</figcaption></figure>

## 3. Attraction to the RJ distribution

We can estimate quantitatively the attraction to the RJ equilibrium by defining a 'distance' to the RJ distribution. For this purpose, we define a quantity that reflects the degree of similarity with the RJ distribution in mode space:

$$D = \frac{\sum_p |n_p^{\text{exp}} - n_p^{\text{RJ}}|}{\sum_p n_p^{\text{exp}} + n_p^{\text{RJ}}}. \qquad (4)$$

Note in particular that this quantity is normalized, $0 \le D \le 1$. We report in Fig. 3(b) the distance $D$ computed for the experimental data $n_p^{\text{exp}}$ at the fiber input (blue), and at the fiber output (red), for different energies $E$. The strong reduction of the distance $D$ from the fiber input to the fiber output clearly confirms the process of thermalization to the RJ equilibrium. It is important to note that this attraction process takes place over a broad range of variation of the energy $E$. Note that, as the energy $E$ increases and approaches the critical value of the transition to condensation $E_{\text{crit}}$, the input-to-output reduction of $D$ decreases as well, a feature that becomes apparent in Fig. 2 (left column), where we can notice that the input and output distributions become similar as the energy increases.

# III. SECOND EXPERIMENT: BEAM CLEANING THROUGH INCIPIENT THERMALIZATION

## A. Differences with the first experiment

There are two differences that distinguish the 1st experimental configuration with the 2nd configuration. (i)
In the 2nd experiment the MMF length is reduced by a factor 4 to $L = 3$m – accordingly the power is increased by the same factor to $N = 28$kW, so as to maintain constant the effective number of nonlinear lengths $L/L_{nl}$ in both experiments. Because of the short fiber length considered, the impact of disorder (weak random coupling) is significantly reduced in this 2nd configuration. (ii) We do not pass the laser beam through a diffuser before injection into the MMF, i.e., we launch a Gaussian-shaped coherent beam (at approximate normal incidence) into the fiber. Notice that the absence of the diffuser in this 2nd configuration does not allow us to make an average over the realizations. The energy $E$ is varied by varying the radius of the launched Gaussian beam: For a given power $N$, by increasing the radius of the Gaussian beam, the population of modes at the fiber input is increased, which in turn increases the energy $E$.

## B. Deviation from the RJ distribution

Following a procedure similar to that carried out in the 1st experiment, we report in Fig. 4 the modal distribution of the power $\tilde{n}_g^{\text{exp}}$, and the energy $\tilde{E}_g^{\text{exp}}$, observed in the 2nd experiment, at the fiber input (blue), and the fiber output (red). We compare the experimental data with the theoretical RJ distribution (green line). We can note that the low-order modes are efficiently attracted toward the RJ equilibrium distribution, while the higher-order modes exhibit some significant deviation to the RJ equilibrium. This aspect becomes even more apparent through the analysis of the distribution of the kinetic energy $\tilde{E}_g^{\text{exp}}$ reported in Fig. 4(c). We can remark that mode groups containing the radially symmetric modes (with $g$ even) are more populated than the other (with $g$ odd), as a natural consequence of the injection of a radially symmetric coherent beam in the MMF. Aside from symmetry considerations, note that nonlinear modal interactions can have different efficiencies also due to the large variation in nonlinear overlap coefficients among the modes, which can significantly affect the rate of convergence to equilibrium [37]. Finally note that the deviation from the RJ distribution is also evidenced by the computation of the distance $D$. As illustrated in Fig. 4(a), the input-to-output decrease of $D$ shows that the system approaches the RJ distribution, although a significant deviation from that distribution is visible at the fiber output.

Although the 2nd experimental configuration only shows an incipient process of RJ thermalization, it evidences a process of beam self-cleaning. This is illustrated in Fig. 5, which reports the modal distribution $\tilde{n}_g^{\text{exp}}$ at low-power, and at high-power, at the output of the MMF. The corresponding 2D intensity distributions at small and high power shows an effect of spatial beam cleaning that is dominated by the behavior of the low order modes, although the higher-order modes significantly deviate from the RJ equilibrium distribution.

<figure><img src="image_6.png" /><figcaption>FIG. 5: **2nd Experiment: Beam self-cleaning.** Modal distribution $\tilde{n}_g$ vs mode group $g$, in the presence of a small power $N = 0.23$kW (linear regime) (a), and high power $N = 28$kW (nonlinear regime) (b). The green line reports the RJ equilibrium distribution. The insets show the 2D intensity distributions: measured (left) and reconstructed (right) output beam profile.</figcaption></figure>

## IV. INTERPRETATION AND CONCLUSION

We have studied the impact of disorder on the process of thermalization in multimode optical fibers by comparing two different experiments realized with a long fiber and an incoherent excitation (1st experiment), and a short fiber with a coherent excitation (2nd experiment). The results of the 1st experiment demonstrate without ambiguity the process of thermalization to the RJ equilibrium distribution, over a broad range of the energy $E$. For small values of $E$, the strong condensate fraction ($n_0 \gg n_{p \neq 0}$) entails a process of beam-cleaning during the propagation through the MMF. In the 2nd experiment, an effect of beam-cleaning can be identified, which can be ascribed to the RJ thermalization of low-order modes, while high-order modes can exhibit a deviation from the RJ distribution.

### A. Discussion of previous results

Previous experimental results revealed a good agreement with the theoretical RJ distribution, especially for low-order modes [42, 44]. On the other hand, some deviation between the experimental measurements and the theoretical RJ distribution was evidenced for high-order modes when the spatial azimuthal complexity grows, see Fig. 8 in Ref. [44]. The results in Ref. [44] were obtained with a short fiber length and the injection of coherent laser beams populating radially symmetric modes, a feature which seems to corroborate the results of our 2nd experiment.

A good agreement has been also obtained between the measured modal distribution and the theoretical RJ distribution for low-order groups of fiber modes [42]. The
<figure><img src="image_7.png" /><figcaption>FIG. 6: **Quasi-resonances accounting for the truncated potential.** Computation of the most favourable non-trivial resonances involving three groups of modes $Δβ_g = \min(2β_g - β_{g-1} - β_{g+1})$, and corresponding resonance efficiency $|Δβ_g|L_{nl}$, where $β_g$ are computed by taking into account the truncation of the parabolic potential due to the fiber cladding. The modes of the fiber can find efficient quasi-resonances $|Δβ_g|L_{nl} \lesssim 1$ that contribute to the thermalization process observed in the 1st experiment, see Fig. 2.</figcaption></figure>

discrepancy between the experiments and the RJ theory for high-order modes was attributed to the lack of efficient resonances among the modes. Indeed, in an ideal parabolic potential, the regular spacing of the mode groups provides efficient phase-matching resonances. In practice however the potential is truncated by the presence of the cladding, which distorts the ladder of mode-group eigenvalues, especially for higher-order modes that overlap more strongly with the cladding. We have computed numerically the mode eigenvalues accounting for the truncation of the parabolic potential for the MMF used in our experiments. We have considered the efficiency of the most favourable non-trivial resonances involving three groups of modes $Δβ_g = \min(2β_g - β_{g-1} - β_{g+1})$, where we recall that *g* indexes the mode group. Then we compute the corresponding resonance efficiency $|Δβ_g|L_{nl}$, where $L_{nl}$ denotes the nonlinear length in the 1st experiment $L_{nl} \approx 0.3$ m. Only quasi-resonances verifying $|Δβ_g|L_{nl} \lesssim 1$ contribute to the thermalization process, while four-mode interactions such that $|Δβ_g|L_{nl} \gg 1$ are non-resonant and do not contribute. We report in Fig. 6 the most favourable non-trivial resonances for each fiber mode, with $g = 1, ..., 13$ (remember $g_{max} = 15$). We can see in Fig. 6 that even higher-order modes of the fiber can find an efficient non-trivial quasi-resonance verifying $|Δβ_g|L_{nl} \lesssim 1$, which indicates that the finite truncation of the parabolic potential should not prevent the thermalization of high-order modes. This corroborates the results reported in the 1st experiment where RJ thermalization is observed over all groups of fiber modes. We also note that the validity of the results reported throughout this paper are not affected by the correction due to the finite truncation of the parabolic potential, i.e., plotting Figs. 1-4 with the correction does not provide visible differences in the figures.

<figure><img src="image_8.png" /><figcaption>FIG. 7: **Numerical simulations:** Simulations of the NLS equation showing the evolutions of the modal components $n_p(z)/N$ (left column), and corresponding evolution of the distance $D(z)$ to the RJ distribution (right column). Fundamental mode $p = 0$ (red), $p = 1$ (dark blue), $p = 2$ (blue), $p = 3$ (light blue), $p = 4$ (cyan). 1st row (a)-(b): Simulation without disorder (weak random mode coupling) and starting from a coherent initial condition. 2nd row (c)-(d): Same simulation as (a)-(b), except that the initial condition is a speckle beam. In (a-d), the system exhibits a coherent regime of mode interaction featured by persistent oscillations that prevent RJ thermalization. 3rd row (e)-(f): Same simulation as (a)-(b), except that disorder is included during the propagation, which corresponds to the 2nd experimental configuration. 4th row (g)-(h): Same simulation as (e)-(f), except that the initial condition is a speckle beam and that the power has been reduced by a factor 4, i.e., 12m of propagation in the bottom row corresponds to 3m in 3rd row – the bottom row corresponds to the 1st experimental configuration. The impact of disorder is more effective in the bottom row (1st experimental configuration) as compared to the 3rd row (2nd experimental configuration): Disorder accelerates RJ thermalization to equilibrium, as revealed by the evolutions of $D(z)$. The horizontal dashed lines (left column) denote the condensate fractions at RJ equilibrium. See the text for parameters.</figcaption></figure>
## B. Interpretation and simulations

The substantial differences that distinguish the results of the 1st and 2nd experiments reported in this work may be interpreted along the lines of the theoretical works in Refs.[35, 37]. We illustrate this aspect in Fig. 7, which reports the numerical simulations of the spatial NLS equation in different configurations. It was shown in Refs.[35, 37] that, by ignoring the impact of disorder (weak random mode coupling due to polarization fluctuations), 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. This regime is illustrated in Fig. 7 (top row), which reports the evolution of the modal components without disorder with a coherent initial condition. Note that a similar regime featured by a phase-sensitive “coherent” modal interaction is also observed by starting the simulation with a speckle beam, see the second row in Fig. 7. The oscillatory behavior of the modal components and the slowing down of RJ thermalization can be related to the existence of Fermi-Pasta-Ulam recurrences, as recently discussed in Ref.[66, 67] in the framework of the weakly nonlinear regime of the 2D NLS equation with a parabolic trapping potential. On the other hand, we recall that a key ingredient for the occurrence of RJ thermalization is the existence of a random phase dynamics among the modes [3, 68] (see also [69, 70]). In other words, a mechanism that breaks the coherent phase-dynamics among the modes is required in order to achieve thermalization in the weakly nonlinear regime considered in the experiments. A natural mechanism is provided by the structural disorder, which is inherent to light propagation in MMFs and whose leading order contribution is associated to random polarization fluctuations [48–51]. This weak random coupling regime is sufficient to introduce a dephasing dynamics among the modes. This is illustrated in the 3rd row of Fig. 7, that reports the same simulation (same coherent initial condition and same parameters) as in the top row, except that a weak random mode coupling has been introduced, with a correlation length $\ell_c = 0.3\text{m}$, and a strength of disorder $2\pi/\sigma = 2.1\text{m}$, see Ref.[35, 37] for details. The simulation in the 3rd row of Fig. 7 then corresponds to the 2nd experimental configuration. We note in Fig. 7 that the oscillatory behavior of the modal components associated to the coherent modal phase dynamics is less pronounced in the 3rd row than in the 1st or 2nd rows: The evolution of the population of the fundamental mode $n_0(z)$ keeps an oscillating behavior, although such a behavior tends to slowly approach the RJ equilibrium value (horizontal dashed black line, $n_0^{\text{RJ}}/N \simeq 0.65$). This observation is corroborated by the evolution of the distance $D(z)$ to the RJ equilibrium that also keeps an oscillatory behavior, though it tends to lower during the propagation (panel (f) of Fig. 7), a feature which is in contrast with

the corresponding evolution without disorder (panels (b), or (d) for $z > 1\text{m}$, in Fig. 7).

This observation is confirmed by the numerical simulation realized in the framework of the 1st experimental configuration, see bottom row in Fig. 7 with an injected speckle beam. In this case a larger fiber length is considered by decreasing the power by a factor 4, so as to keep constant the effective number of nonlinear interaction lengths, i.e., 12m of propagation in the bottom row of Fig. 7 corresponds to 3m in the 3rd row (or in the 1st and 2nd rows). Because of the larger fiber length, the impact of disorder is more effective for the 1st experimental configuration (bottom row) as compared to the 2nd configuration (3rd row): The condensate fraction in Fig. 7(g) reaches the RJ thermal value, while the distance to RJ equilibrium tends to decrease to zero in an almost monotonic fashion, see Fig. 7(h). Interestingly, the evolutions of $n_0(z)$ and $D(z)$ in the 1st experimental configuration (bottom row) do not exhibit the oscillatory behaviour evidenced in the 2nd configuration (3rd row), which confirms the absence of modal phase-correlations and the significant acceleration of RJ thermalization due to the disorder.

The simulations reported in Fig. 7 qualitatively reproduce the experimental results. However, it should be stressed that a power larger by a factor three has been considered to accelerate the dynamics in both the 1st and 2nd configurations. Then although the purely spatial model considered in the simulations captures the essential features of the experiments, a quantitative description of the experimental results would require a detailed analysis of the temporal averaging effect inherent to the pulsed regime considered in the experimental measurements [71–73].

We finally note that, on the basis of the wave turbulence theory, we have derived a kinetic equation describing the evolutions of the modal components $n_p(z)$, which revealed that polarization fluctuations significantly accelerate the process of RJ thermalization, see Refs.[35, 37]. This theoretical approach then explains, at a pure qualitative level, the fast process of RJ thermalization observed in the 1st experiment reported in this work. We recall that, at variance with the 2nd experiment, the 1st experiment is characterized by two factors: (i) a larger fiber length that enhances the impact of disorder, (ii) the injection of a speckle beam owing to the diffuser. The combination of these two factors leads to a random phase dynamics among the modes, which favours the process of RJ thermalization observed in the 1st experiment. Finally, this work contributes to the understanding of the fast process of optical thermalization observed in multi-mode optical fibers, a feature that may also be important for possible future extensions to the thermalization of multiple beams [74, 75]. From a broader perspective, this work also contributes to the understanding of the interplay of nonlinearity and disorder [76–80], in relation with the paradigm of statistical light-mode dynamics, glassy behaviors, and complexity science [81].
## V. ACKNOWLEDGMENTS

The authors are grateful to S. Rica, I. Carusotto and V. Doya for fruitful discussions. Fundings: Centre national de la recherche scientifique (CNRS), Conseil régional de Bourgogne Franche-Comté, iXCore Research Fondation, Agence Nationale de la Recherche (ANR-19-CE46-0007, ANR-15-IDEX-0003, ANR-21-ESRE-0040). Calculations were performed using HPC resources from DNUM CCUB (Centre de Calcul, Université de Bourgogne).

## VI. APPENDIX: CONVERGENCE TO RJ EQUILIBRIUM WITH THE NUMBER OF SPECKLE REALIZATIONS

In this Appendix, we show that the modal distribution recorded experimentally $n_p^{\text{exp}} = \langle |a_p^{\text{exp}}|^2 \rangle$ converges toward the expected RJ distribution $n_p^{\text{RJ}}$ as the number of realizations of speckles beams involved in the average is increased. Here, we analyze the results retrieved from the 1st experimental configuration, see Figs. 1-2. The convergence of $n_p^{\text{exp}}$ to the RJ distribution follows qualitatively the behavior expected theoretically.

As discussed in section II B, we have recorded in the experiments 2×1000 realizations of the NF and FF intensity distributions for the same power $N = 7$ kW and different values of the energy $E$. For each individual realization, we have retrieved the modal distribution $|a_p^{\text{exp}}|^2$. We have partitioned such an ensemble of realizations within small intervals of energies $\Delta E/E_0 = 0.125$, with $E_0 = N\beta_0$ is the minimum of the energy. Next, we have performed an average over the realizations of the modal distributions for each specific energy interval $\Delta E$. In the following, we consider a particular energy interval around $E/(N\beta_0) = 6$ and we study the convergence to the RJ equilibrium distribution with the number of realizations of speckle beams. We define the distance between the experimental data and the RJ equilibrium as follows:

$$ \mathcal{D}_Q = \frac{\sum_p |n_p^{\text{exp},Q} - n_p^{\text{RJ}}|}{\sum_p n_p^{\text{exp},Q} + n_p^{\text{RJ}}}. \qquad (5) $$

where the experimental distribution $n_p^{\text{exp},Q} = (1/Q) \sum_{j=1}^Q |a_p^{\text{exp},j}|^2$ is the empirical average over $Q$ realizations of the modal distributions $|a_p^{\text{exp},j}|^2$. We report the distance $\mathcal{D}_Q$ vs number of realizations $Q$ in Fig. 8(a).

We compare the experimental results with the theoretical expression of the corresponding distance that we

have derived in Ref. [55] (Eq. (3) in the Supplementary Methods):

$$ \mathcal{D}_Q \approx \frac{Q^{Q-1}}{(Q-1)!} e^{-Q}. \qquad (6) $$

<figure><img src="image_9.png" /><figcaption>FIG. 8: Experimental convergence to the RJ distribution with the number of realizations. (a) The blue points denote the distance $\mathcal{D}_Q$ [Eq. (5)] between the experimental modal distribution $n_p^{\text{exp},Q}$ averaged over $Q$ realizations of speckle beams, and the theoretical RJ equilibrium distribution $n_p^{\text{RJ}}$. The black line reports the theoretical distance $\mathcal{D}_Q$ computed analytically, from Eq. (6). Note that $\mathcal{D}_Q$ is normalized, $0 \le \mathcal{D}_Q \le 1$. (b) Experimental modal distribution $n_p^{\text{exp},Q}$ averaged over $Q$ realizations [$Q = 1$ (blue), $Q = 10$ (green), $Q = 50$ (red)]. Corresponding theoretical RJ distribution $n_p^{\text{RJ}}$ (black line).</figcaption></figure>

For $Q \ge 8$, we have $\mathcal{D}_Q \approx 1/\sqrt{2\pi Q}$. Note that Eq. (6) has been found in quantitative agreement with the numerical simulations (see Fig. 4 in the Supplementary Material of Ref. [55]). We can note that the distance $\mathcal{D}_Q$ to the RJ distribution computed from the experimental results (blue points) decreases with the number of realizations by following qualitatively the theoretical behavior (solid black line). The theoretical curve given by Eq. (6) assumes that the Gerchberg-Saxton algorithm, as well as the sampling of the cameras used in the experiments do not introduce errors. The impact of the two errors can explain that the distance $\mathcal{D}_Q$ for the experimental data reported in Fig. 8(a) is slightly above the theoretical prediction. Finally, we report in Fig. 8(b) the experimental modal distributions $n_p^{\text{exp},Q}$ averaged over $Q$ realizations of speckle beams, which shows the convergence toward the expected RJ distribution $n_p^{\text{RJ}}$ as $Q$ increases.

[1] R.Z. Sagdeev, D.A. Usikov, and G.M. Zaslavsky, *Non-linear Physics* (Harwood Publishing, Chur, Switzerland, 1988).

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

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

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

[5] S. K. Turitsyn, S. A. Babin, E. G. Turitsyna, G. E. Falkovich, E. Podivilov, D. Churkin, Optical Wave Turbulence, in Advances in Wave Turbulence (World Scientific, Singapore, 2013).

[6] D. Churkin, I. Kolokolov, E. Podivilov, I. Vatnik, S. Vergeles, I. Terekhov, V. Lebedev, G. Falkovich, M. Nikulin, S. Babin, S. Turitsyn, Wave kinetics of a random fibre laser, Nature Communications **2**, 6214 (2015).

[7] V. Zakharov, A. Pushkarev, V. Shvets, V. Yan'kov, Soliton turbulence, JETP Lett. **48**, 83 (1988).

[8] R. Jordan, C. Josserand, Self-organization in nonlinear wave turbulence, Phys. Rev. E **61**, 1527 (2000).

[9] B. Rumpf, A.C. Newell, Coherent structures and entropy in constrained, modulationally unstable, nonintegrable systems, Phys. Rev. Lett. **87**, 054102 (2001).

[10] B. Rumpf, A.C. Newell, Localization and coherence in nonintegrable systems, Physica D **184**, 162 (2003).

[11] V.E. Zakharov, F. Dias, A. Pushkarev, One-dimensional wave turbulence, Phys. Reports **398**, 1 (2004).

[12] J. Laurie, U. Bortolozzo, S. Nazarenko, S. Residori, One-Dimensional Optical Wave Turbulence: Experiment and Theory, Phys. Reports **514**, 121 (2012).

[13] B. Rumpf, A.C. Newell, V.E. Zakharov, Turbulent Transfer of Energy by Radiating Pulses, Phys. Rev. Lett. **103**, 074502 (2009).

[14] M. Conforti, A. Mussot, J. Fatome, A. Picozzi, S. Pitois, C. Finot, M. Haelterman, B. Kibler, C. Michel, and G. Millot, Turbulent dynamics of an incoherently pumped passive optical fiber cavity: Quasisolitons, dispersive waves, and extreme events, Phys. Rev. A **91**, 023823 (2015).

[15] S. Pickartz, U. Bandelow, S. Amiranashvili, Adiabatic theory of solitons fed by dispersive waves, Phys. Rev. A **94**, 033811 (2016).

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

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

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

[19] S. Nazarenko, M. Onorato, Wave turbulence and vortices in Bose-Einstein condensation, Physica D **219**, 1 (2006).

[20] G. Düring, A. Picozzi, S. Rica, Breakdown of weak-turbulence and nonlinear wave condensation, Physica D **238**, 1524 (2009).

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

[22] A. Picozzi, S. Rica, Condensation of classical optical waves beyond the cubic nonlinear Schrödinger equation, Optics Communications **285**, 5440 (2012).

[23] S. Nazarenko, M. Onorato, D. Proment, Bose-Einstein condensation and Berezinskii-Kosterlitz-Thouless transition in the two-dimensional nonlinear Schrödinger model, Phys. Rev. A **90**, 013624 (2014).

[24] E. Turitsyna, G. Falkovich, A. El-Taher, X. Shu, P. Harper, and S. Turitsyn, Optical turbulence and spectral condensate in long fibre lasers, Proc. R. Soc. London Ser. A **468**, 2145 (2012).

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

[26] C. Sun, S. Jia, C. Barsi, S. Rica, A. Picozzi, J.W. Fleischer, Observation of the kinetic condensation of classical waves, Nat. Phys. **8**, 470 (2012).

[27] J. Bloch, I. Carusotto, M. Wouters, Nonequilibrium Bose-Einstein condensation in photonic systems, Nature Physics Reviews **4**, 470 (2022).

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

[29] P. Aschieri, J. Garnier, C. Michel, V. Doya, A. Picozzi, Condensation and thermalization of classsical optical waves in a waveguide, Phys. Rev. A **83**, 033838 (2011).

[30] In the absence of a waveguide, RJ thermalization is not well defined. Consider an experiment realized in a 2D bulk medium, where the RJ distribution takes the form $n_k^{RJ} = T/(\omega_k - \mu)$, with $\omega_k \sim k^2$ the dispersion relation of the NLS equation. The power $N = \int n_k^{RJ} d^2k$ and the kinetic energy $E = \int \omega_k n_k^{RJ} d^2k$ exhibit an ultraviolet divergence at thermal equilibrium [26]. To see this, consider the behaviour of the RJ spectrum for large frequencies $n_k^{RJ} \sim 1/k^2$, then $N \sim \int dk/k$, and $E \sim \int dk$, which are divergent integrals for $k \to \infty$. The finite number of modes of the waveguide introduces an effective frequency cut-off that regularizes such an ultraviolet divergence [29, 40].

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

[32] Z. Liu, L.G. Wright, D.N. Christodoulides, F.W. Wise, Kerr self-cleaning of femtosecond-pulsed beams in graded-index multimode fiber, Optics Letters **41** 3675 (2016).

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

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

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

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

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

[38] F.O. Wu, A.U. Hassan, D.N. Christodoulides, Thermo-dynamic theory of highly multimoded nonlinear optical systems, Nature Photon. **13**, 776 (2019).

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

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

[41] M. Fabert, M. Sapantan, K. Krupa, A. Tonello, Y. Lev entoux, S. Février, T. Mansuryan, A. Niang, B. Wetzel, G. Millot, S. Wabnitz, V. Couderc, Coherent combining of self-cleaned multimode beams, Scient. Rep. **10**, 20481 (2020).

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

[43] K. Baudin, A. Fusaro, J. Garnier, N. Berti, K. Krupa, I. Carusotto, S. Rica, G. Millot, A. Picozzi, Energy and wave-action flows underlying Rayleigh-Jeans thermalization of optical waves propagating in a multimode fiber, Europhys. Lett. **134**, 14001 (2021).

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

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

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

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

[48] A. Mecozzi, C. Antonelli, M. Shtaif, Nonlinear propagation in multimode fibers in the strong coupling regime, Opt. Exp. **20**, 11673 (2012).

[49] A. Mecozzi, C. Antonelli, M. Shtaif, Coupled Manakov equations in multimode fibers with strongly coupled groups of modes, Opt. Exp. **20**, 23436 (2012).

[50] S. Mumtaz, R.J. Essiambre, G.P. Agrawal, Nonlinear propagation in multimode and multicore fibers: Generalization of the Manakov equations, J. Lightw. Technol. **31**, 398 (2013).

[51] Y. Xiao, R.-J. Essiambre, M. Desgroseilliers, A.M. Tulino, R. Ryf, S. Mumtaz, G.P. Agrawal, Theory of intermodal four-wave mixing with random linear mode coupling in few-mode fibers, Opt. Exp. **22**, 32039 (2014).

[52] O.S. Sidelnikov, E.V. Podivilov, M.P. Fedoruk, S. Wabnitz, Random mode coupling assists Kerr beam self-cleaning in a graded-index multimode optical fiber, Optical Fiber Technology, **53**, 101994 (2019).

[53] D. Psaltis, C. Moser, Imaging with multimode fibers, Opt. and Photon. News **27**, 24 (2016).

[54] P. Caramazza, O. Moran, R. Murray-Smith, D. Faccio,

Transmission of natural scene images through a multimode fibre, Nat Commun **10**, 2029 (2019).

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

[56] E.S. Manuylovich, V.V. Dvoyrin, S.K. Turitsyn, Fast mode decomposition in few-mode fibers, Nature Communications **11**, 5507 (2020).

[57] O. Shapira, A.F. Abouraddy, J.D. Joannopoulos, Y. Fink, Complete modal decomposition for optical waveguides, Phys. Rev. Lett. **94**, 143902 (2005).

[58] Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, Phase retrieval with application to optical imaging: A contemporary overview, IEEE Signal Process. Mag. **32**, 87 (2015).

[59] J.R. Fienup, Reconstruction of an object from the modulus of its Fourier transform, Opt. Lett. **3**, 27 (1978).

[60] J.R. Fienup, Phase retrieval algorithms: a comparison, Applied Optics **21**, 2758 (1982).

[61] S. Pitois, S. Lagrange, H.R. Jauslin, A. Picozzi, Velocity locking of incoherent nonlinear wave packets, Phys. Rev. Lett. **97**, 033902 (2006).

[62] S. Lagrange, H.R. Jauslin, A. Picozzi, Thermalization of the dispersive three-wave interaction, Europhys. Lett. **79**, 64001 (2007).

[63] B. Barviau, B. Kibler, A. Kudlinski, A. Mussot, G. Millot, A. Picozzi, Experimental signature of optical wave thermalization through supercontinuum generation in photonic crystal fiber, Optics Express **17**, 7392 (2009).

[64] B. Barviau, B. Kibler, A. Picozzi, Wave-turbulence approach of supercontinuum generation: Influence of self-steepening and higher-order dispersion Phys. Rev. A **79**, 063840 (2009).

[65] F.O. Wu, Q. Zhong, H. Ren, P.S. Jung, K.G. Makris, D.N. Christodoulides, Thermalization of Light's Orbital Angular Momentum in Nonlinear Multimode Waveguide Systems, Phys. Rev. Lett. **128**, 123901 (2022).

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

[67] A. Biasi, O. Evnin, B.A. Malomed, Obstruction to ergodicity in nonlinear Schrödinger equations with resonant potentials, arXiv:2304.10308 (2023).

[68] S. Chibbaro, G. Dematteis, C. Josserand, L. Rondoni, Wave-turbulence theory of four-wave nonlinear interactions, Phys. Rev. E **96**, 021101(R) (2017).

[69] M. Guasoni, J. Garnier, B. Rumpf, D. Sugny, J. Fatome, F. Amrani, G. Millot, A. Picozzi, Incoherent Fermi-Pasta-Ulam Recurrences and Unconstrained Thermalization Mediated by Strong Phase Correlations, Phys. Rev. X **7**, 011025 (2017).

[70] J. Zaleski, M. Onorato, and Y.V. Lvov, Anomalous Correlators in Nonlinear Dispersive Wave Systems, Phys. Rev. X **10**, 021043 (2020).

[71] Y. Leventoux, G. Granger, K. Krupa, A. Tonello, G. Millot, M. Ferraro, F. Mangini, M. Zitelli, S. Wabnitz, S. Février, V. Couderc, 3D time-domain beam mapping for studying nonlinear dynamics in multimode optical fibers, Optics Letters **46**, 66 (2021).

[72] Y. Leventoux, G. Granger, K. Krupa, T. Mansuryan, M. Fabert, A. Tonello, S. Wabnitz, V. Couderc, S. Février, Frequency-resolved spatial beam mapping in multimode
fibers: application to mid-infrared supercontinuum gen-
eration, Optics Letters **46**, 3717 (2021).

[73] M. Conforti, C. Mas Arabi, A. Mussot, and A. Kudlinski, Fast and accurate modeling of nonlinear pulse propagation in graded-index multimode fibers, Optics Letters **42**, 4004 (2017).

[74] M. Ferraro, F. Mangini, Y. Leventoux, A. Tonello, M. Zitelli, T. Mansuryan, Y. Sun, S. Fevrier, K. Krupa, D. Kharenko, S. Wabnitz, V. Couderc, Multimode optical fiber beam-by-beam cleanup, Journal of Lightwave Technology, 1 (2023).

[75] M. Ferraro, F. Mangini, F.O. Wu, M. Zitelli, D.N. Christodoulides, S. Wabnitz, Calorimetry of photon gases in nonlinear multimode optical fibers, arXiv:2212.12781

[76] N. Cherroret, T. Karpiuk, B. Grémaud, C. Miniatura, Thermalization of matter waves in speckle potentials, Phys. Rev. A **92**, 063614 (2015).

[77] T. Scoquart, P.-E. Larré, D. Delande, N. Cherroret, Weakly interacting disordered Bose gases out of equilibrium: From multiple scattering to superfluidity, European Physics Letters **132**, 66001 (2020).

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

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

[80] P. Haldar, S. Mu, B. Georgeot, J. Gong, C. Miniatura, G. Lemarié, Prethermalization and wave condensation in a nonlinear disordered Floquet system, arXiv:2109.14347

[81] C. Conti, E. DelRe, Photonics and the Nobel Prize in Physics, Nature Photonics **16**, 6 (2022).