# Observation of light thermalization to negative-temperature Rayleigh-Jeans equilibrium states in multimode optical fibers

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

► To cite this version:

K. Baudin, J. Garnier, A. Fusaro, N. Berti, C. Michel, et al., Observation of light thermalization to negative-temperature Rayleigh-Jeans equilibrium states in multimode optical fibers. Physical Review Letters, 2023, 130 (6), pp.063801. 10.1103/PhysRevLett.130.063801. hal-04198191

HAL Id: hal-04198191

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

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>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Observation of light thermalization to negative temperature Rayleigh-Jeans equilibrium states in multimode optical fibers

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

¹ Laboratoire Interdisciplinaire Carnot de Bourgogne, CNRS, Université Bourgogne Franche-Comté, 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

Although the temperature of a thermodynamic system is usually believed to be a positive quantity, under particular conditions, negative temperature equilibrium states are also possible. Negative temperature equilibriums have been observed with spin systems, cold atoms in optical lattices and two-dimensional quantum superfluids. Here we report the observation of Rayleigh-Jeans thermalization of light waves to negative temperature equilibrium states. The optical wave relaxes to the equilibrium state through its propagation in a multimode optical fiber, i.e., in a conservative Hamiltonian system. The bounded energy spectrum of the optical fiber enables negative temperature equilibriums with high energy levels (high order fiber modes) more populated than low energy levels (low order modes). Our experiments show that negative temperature speckle beams are featured, in average, by a non-monotonous radial intensity profile. The experimental results are in quantitative agreement with the Rayleigh-Jeans theory without free parameters. Bringing negative temperatures to the field of optics opens the door to the investigation of fundamental issues of negative temperature states in a flexible experimental environment.

**Introduction.**- Temperature is a central concept of statistical mechanics and often reflects a measure of the amount of disordered motion in a classical ideal gas. Although this intuitive notion is correct for many physical systems, one should keep in mind that the concept of temperature is by far more subtle. A detailed analysis of the concept of temperature, and of its relationship with energy and entropy shows that, under suitable conditions, the entropy can decrease with the energy, thus allowing for the existence of equilibrium states at negative temperatures (NT). Starting from the seminal works by Onsager [1] and Ramsey [2], who originally conceived the physical idea and the first theoretical approaches, during the last decades, many works have been devoted to the theoretical understanding of these unusual equilibrium states. Despite the fact that the existence of a NT equilibrium has created its own share of confusion in relation with the definition of the entropy [3, 4], NTs are now broadly accepted in line with different experimental observations [5–12]. NTs were originally observed experimentally in nuclear spin systems [13]. More recently, NTs were observed with cold atoms in optical lattices [14]. Furthermore, NTs originally predicted by Onsager in the statistical description of point vortices [1] have been recently observed in 2D quantum superfluids [15, 16].

In this Letter we present an experimental optical setup in which we report the observation of light thermalization to NT equilibrium states. Our system is based on the nonlinear propagation of speckle beams in a multimode optical fiber (MMF). Because of the presence of a finite number of modes supported by the MMF, the spectrum exhibits both lower and upper bounds for the

energy levels. The bounded spectrum, combined to the nonlinear four-wave interaction, are responsible for the process of Rayleigh-Jeans thermalization to NT equilibrium states [17, 18]. We stress that, at variance with other experiments where photon thermalization is driven by a thermal heat bath [19–22], here light thermalization takes place in a conservative Hamiltonian system. RJ thermalization to usual *positive temperature* equilibriums has been recently demonstrated experimentally in MMFs [23–26], on the basis of a spatial beam-cleaning effect [27–31]. As described by the wave turbulence theory [32–36] applied to MMFs [37–40], the thermalization to a positive temperature equilibrium is characterized by a transfer of power (particle number) toward the low-order modes of the MMF. In marked contrast, here we report the observation of thermalization to a NT equilibrium featured by a power transfer to high-order modes (direct flow of particles), as well as a transfer of energy to low-order modes (inverse flow of energy). Consequently, the NT equilibrium is characterized by an inverted modal population, in which high-order modes are more populated than low-order modes.

Our experimental optical setup can be used as a simple and flexible testbed to explore fundamental issues related to NT states that are discussed in conclusion, e.g., Carnot cycles operating between temperatures of opposite signs, or inverted turbulence cascades featured by an analogue process of condensation at NT.

**Experimental system.**- The experiment is based on the single pass propagation of speckle beams through a MMF. The subnanosecond pulses delivered by a Nd:YAG laser ($\lambda = 1.06\mu$m) are transmitted through a spiral
phase plate and then through a diffuser before injection of the speckle beam into a 12m long graded-index MMF (i.e., parabolic-shaped trapping potential), which guides $M = 45$ modes, i.e., nine groups of degenerate modes. The energy levels (fiber eigenvalues) are well approximated by the ones of an harmonic potential $\beta_p = \beta_0(p_x + p_y + 1)$, where $\{p\}$ labels the two integers $(p_x, p_y)$ that specify a mode (see Supplementary Material). We denote by $|a_p|^2$ the power in the mode $p$, with the total power $N = \sum_p |a_p|^2$ [44].

The experiment is realized in the weakly nonlinear regime, where linear effects dominate over nonlinear effects $L_{lin} \sim \beta_0^{-1} \sim 0.1 \text{mm} \ll L_{nl} = 1/(\gamma N) \sim 20 \text{cm}$, $\gamma$ being the nonlinear coefficient of the MMF. Accordingly, we do not consider NT states associated to nonlinear coherent structures, e.g., breathers [10, 45, 46]. Since $L_{lin} \ll L_{nl}$, we only retain the linear contribution to the Hamiltonian, $E = \sum_p \beta_p |a_p|^2$ [23–25]. We have verified the conservation of the power $N$ and the energy $E$ through propagation in the NT region for each realization of a speckle beam, which confirms that the coupling between guided modes and leaky modes of the fiber can be neglected (see Supplementary Material).

RJ thermalization is driven by the four-wave nonlinear interaction through the propagation in the MMF. The speckle beam is expected to relax toward the thermodynamic equilibrium state described by the RJ distribution [17, 23–25, 38, 39]:

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

where $T$ and $\mu$ are the temperature and chemical potential, while $n_p = \langle |a_p|^2 \rangle$ denotes the modal power averaged over the realizations of the speckle beams. We have at equilibrium $N = T \sum_p (\beta_p - \mu)^{-1}$ and $E = T \sum_p \beta_p / (\beta_p - \mu)$, with $(T, \mu)$ uniquely determined by $(N, \bar{E})$ – we deal with a microcanonic description ($T$ is not defined by a thermostat, it is in units of $\text{W} \cdot \text{m}^{-1}$) [36]. Note that the RJ distribution refers to the classical, low-energy, limit of the Bose-Einstein distribution [32], describing highly occupied fiber modes.

**Negative temperatures.**- The irreversible process of RJ thermalization is described by the wave turbulence theory [32–36], which provides a nonequilibrium description of light propagation in MMFs [37–40]. An equilibrium thermodynamic formulation of multimode optical systems has been recently developed [17, 47]. We report in Fig. 1 the relative entropy $S = \sum_p \log(n_p^{\text{RJ}})$ as a function of the energy for the MMF used in our experiments with $g = 9$ groups of degenerate modes. Because the spectrum of the fiber is bounded, $\beta_0 \le \beta_p \le \beta_{\text{max}} = g\beta_0$, the system possesses both lower and upper energy bounds

$$E_{\text{min}} = N\beta_0 \le E \le E_{\text{max}} = N\beta_{\text{max}}. \quad (2)$$

Starting at minimum energy $E_{\text{min}}$, where only the fundamental mode is populated, an increase in energy leads

<figure><img src="image_2.png" /><figcaption>FIG. 1: Negative temperatures and inverted modal population. (a) RJ equilibrium distribution $n_p^{\text{RJ}}$ for positive temperature $T > 0$ ($E < E_*$) where low-order modes are more populated, and negative temperatures $T < 0$ ($E > E_*$) featured by an inverted modal population; while for $1/T \to 0$ ($E = E_*$), $n_p^{\text{RJ}} = \text{const}$. (b) Relative entropy $S$ vs energy $E$, showing that $1/T = (\partial S / \partial E)_{N,M} < 0$ requires $E > E_*$. (c) Temperature $T$ vs energy $E$. Negative temperatures $T < 0$ occur for $E > E_*$ with $E_*/E_{\text{min}} \simeq 6.33$ (vertical dashed black line). The vertical purple lines in (b-c) denote the six values of $E$ considered in Fig. 2.</figcaption></figure>

to an occupation of a larger number of fiber modes and therefore an increase in entropy. As the temperature approaches infinity, all fiber modes become equally populated $n_p^{\text{RJ}} = \text{const}$, and the entropy reaches a maximum for $E = E_* = N \langle \beta_p \rangle = E_{\text{min}}(2g+1)/3$. NT equilibrium states arise for $E > E_*$, where the entropy decreases by increasing the energy, $1/T = (\partial S / \partial E)_{M,N} < 0$. The condition $E > E_*$ can be achieved if high-order modes are more populated than low-order modes. Note that NT equilibrium states persist in the thermodynamic limit (see Supplementary Material).

**RJ thermalization to NT equilibriums.**- At variance with usual experiments of spatial beam cleaning and RJ thermalization [24–30], here we study the thermalization for different values of the energy $E$, while keeping constant the power $N$. Indeed, by passing the laser beam through a diffuser before injection into the fiber, we can vary the amount of randomness of the speckle beam by keeping $N = \text{const}$ – the larger the randomness of the speckle beam, the higher the energy $E$. Accordingly, we study RJ thermalization over a broad range of variation of the energy. In order to further increase the energy beyond the threshold for NT ($E > E_*$), we pass the beam through a spiral phase plate before the diffuser, i.e., we generate a speckle beam from a doughnut-like intensity
<figure><img src="image_3.png" /><figcaption>FIG. 2: Rayleigh-Jeans thermalization to NT equilibriums. Experimental modal distributions averaged over realizations $n_p^{\text{exp}} = \langle |a_p^{\text{exp}}|^2 \rangle$, at the fiber input (blue), at the fiber output (red). Corresponding RJ equilibrium distribution $n_p^{\text{RJ}}$ (green). Note the quantitative agreement between $n_p^{\text{RJ}}$ and the experimental output distribution $n_p^{\text{exp}}$ (red). The six panels correspond to different values of $E$, or equivalently different $T$, see the six vertical purple lines in Fig. 1(b-c). The modal distribution peaked on the lowest mode for $T > 0$ (a), gets inverted for $T < 0$ (b-f). An average over ≈35 realizations of speckle beams is considered for each panel. The fiber modes are sorted from the fundamental one ($\beta_0$) to the highest mode group (nine-fold degenerate with $\beta_{\text{max}} = 9\beta_0$). Degenerate modes are equally populated at equilibrium, leading to a staircase distribution $n_p^{\text{RJ}}$.</figcaption></figure>

distribution, which enables the excitation of higher order fiber modes.

The accurate measurements of the near-field and far-field intensity distributions allowed us to retrieve the modal power distribution $n_p^{\text{exp}}$. To obtain the mode decomposition, several interferometric approaches based on use of a reference beam have been exploited to study light thermalization in MMFs [24–26]. Here, in contrast to the previous works [23–26], we use a non-interferometric numerical mode decomposition procedure that is based on the Gerchberg-Saxton algorithm. It allows us to retrieve the transverse phase profile of the speckle field from the near-field and far-field intensity distributions measured in the experiments [48–51]. By projecting the retrieved complex field over the fiber modes, we get the complete modal distribution.

The RJ distribution being in essence a statistical distribution, its comparison with the experiments requires an average over realizations of speckle beams. We have recorded 2×300 realizations of the near-field and far-field intensity distributions for the same power ($N = 7\,\text{kW}$)

<figure><img src="image_4.png" /><figcaption>FIG. 3: Energy flows in mode space. Experimental energy distributions averaged over 50 realizations $\epsilon_p^{\text{exp}} = \beta_p n_p^{\text{exp}}$, at the fiber input (blue), output (red). The arrow indicates the energy flow to low-order modes. Corresponding RJ equilibrium distribution $\epsilon_p^{\text{RJ}} = \beta_p n_p^{\text{RJ}}$ (green line), which is in quantitative agreement with the experimental output distribution (red).</figcaption></figure>

and different energies $E$. For each individual speckle realization, we retrieve the modal distribution $|a_p^{\text{exp}}|^2$. We partition the ensemble of 300 realizations of $\{|a_p^{\text{exp}}|^2\}$ within small energy intervals $[E - \Delta E, E + \Delta E]$ with $\Delta E = 0.125E_{\text{min}}$. We perform an average over the realizations of the modal distributions for each energy interval, which provides the averaged modal distribution $n_p^{\text{exp}} = \langle |a_p^{\text{exp}}|^2 \rangle$. This procedure is applied at the fiber output ($L = 12\,\text{m}$), and fiber input (after 20cm of propagation). The error in the procedure has been computed theoretically and numerically, it decreases with the number of realizations and has been found remarkably small (relative standard deviations of ≈ 6%), see Supplementary Material.

We report in Fig. 2 the averaged modal distributions $n_p^{\text{exp}}$ at the fiber input (blue) and output (red), for different values of the energy $E$, or equivalently the temperature $T$ (purple lines in Fig. 1(b-c)). The data are compared with the theoretical RJ distribution $n_p^{\text{RJ}}$. We stress that there are no adjustable parameters between $n_p^{\text{exp}}$ and $n_p^{\text{RJ}}$: The parameters ($T, \mu$) in $n_p^{\text{RJ}}$ are uniquely determined by $N$ and $E$ measured in the experiments. We observe in Fig. 2 an excellent agreement between $n_p^{\text{exp}}$ (red circles) and $n_p^{\text{RJ}}$ (green line), for both $T > 0$ and $T < 0$. Fig. 2 then shows that NT equilibriums constitute *attractor states* for the random wave, whose robustness has a thermodynamic origin – maximum entropy state for a given pair $(N, E)$.

*Energy flows in mode space.*- The conventional thermalization to positive temperatures is characterized by an energy flow to high-order modes [33, 34, 38]. Thermalization to NTs typically occurs through an inverse energy flow to low-order modes [18]. This is illustrated in Fig. 3, which shows that the energy distribution $\epsilon_p^{\text{exp}} = \beta_p n_p^{\text{exp}}$ at low-order modes increases through propagation in the MMF and reaches the theoretical RJ equilibrium distribution $\epsilon_p^{\text{RJ}} = \beta_p n_p^{\text{RJ}}$.
<figure><img src="image_5.png" /><figcaption>FIG. 4: Oscillating radial intensity distribution at NT. Averaged intensity distribution $I^{\text{exp}}(|\mathbf{r}|)$ as a function of the radial (angle-averaged) distance $|\mathbf{r}|$ (red). Note the quantitative agreement with the theoretical RJ intensity distribution $I^{\text{RJ}}(|\mathbf{r}|)$ in Eq.(3) (dashed green). The oscillating behavior of the intensity distribution is a signature of the NT equilibrium. Inset: corresponding 2D intensity averaged over the realizations (the radius of the circle is the fiber radius).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5: Macroscopic population of the highest energy level. Fraction of power $\tilde{n}_g/N$ into the highest mode group $g = 9$ vs energy $E/E_{\text{min}}$. Experimental measurements at the fiber output: The blue circles refer to the linear regime (small power), the red circles to the nonlinear regime (high power). The green line denotes the RJ equilibrium theory. By increasing the energy, the power goes to the highest energy level, $\tilde{n}_g/N \to 1$ as $E/E_{\text{min}} \to 9$.</figcaption></figure>

**Oscillating radial intensity distribution.**- The intensity distribution $I^{\text{RJ}}(|\mathbf{r}|)$ of usual positive temperature equilibriums is, in average, a monotonic decreasing function with the radial distance $|\mathbf{r}|$ [23]. This is consistent with the intuitive idea that low-order modes localized near-by the fiber center are the most populated ones. In marked contrast, the inverted modal population of NT equilibriums are characterized by an oscillating behaviour of the radial intensity distribution. This is illustrated in Fig. 4, which reports the averaged radial intensity distribution $I^{\text{exp}}(|\mathbf{r}|)$ (with $\Delta E = 0.25 E_{\text{min}}$, $E/E_{\text{min}} = 7.9$). The theoretical RJ intensity distribution reads

$$I^{\text{RJ}}(\mathbf{r}) = \sum_{p} n_{p}^{\text{RJ}} u_{p}^{2}(\mathbf{r}), \quad (3)$$

where $u_p(\mathbf{r})$ denotes the fiber modes [23]. The number of radial oscillations in Fig. 4 is given by the most oscillating mode of the fiber, namely the mode LP$_{04}$ that exhibits 5 oscillations.

**Experiments by increasing power.**- We have studied the optical field at the output of the MMF, with a small power $N = 0.23\,\text{kW}$ (linear regime), and a high power $N = 7\,\text{kW}$ (nonlinear regime). Since the MMF length is kept fixed ($L = 12\,\text{m}$), the effective number of nonlinear interaction lengths increases by increasing the power. Fig. 5 reports the fraction of power that populates the highest group of degenerate modes of the MMF, $\tilde{n}_g/N$ for $g = 9$. The output field (red) reaches the equilibrium RJ theory (green line) in the nonlinear regime. The highest energy level gets macroscopically populated by increasing the energy, or equivalently by increasing the temperature of negative sign (see Fig. 1).

**Conclusion and perspectives.**- We have reported the observation of RJ thermalization to NT equilibrium states through light propagation in graded-index MMFs. This non-equilibrium process of NT thermalization can be described by a wave turbulence kinetic equation, which is

found in agreement with the simulations of the nonlinear Schrödinger equation (see Supplementary Material). Our NT experiment then paves the way for the study of Zakharov-Kolmogorov turbulence cascades [32–34] that are inverted with respect to those underlying usual positive temperature thermalization (e.g., inverse energy flow in Fig. 3).

Along this line, our work suggests a previously unrecognized process of inverted condensation at NTs: At variance with usual condensation at positive temperature where the lowest energy level gets macroscopically populated by decreasing the temperature ($T \to 0^+$, or $E \to E_{\text{min}}$) [23, 33–36], at NT an inverted condensation process occurs into the highest energy level as the temperature increases to zero ($T \to 0^-$, or $E \to E_{\text{max}}$). While we provide a preliminary study of this effect through the macroscopic population of the highest energy level (Fig. 5), the observation of the transition to condensation requires MMFs with larger number of modes (see Supplementary Material).

In our work NT states are obtained directly, which is in contrast with magnetic systems and cold atoms where the excitation of NT states requires first the creation of a positive temperature state and then its subsequent inversion through suitable procedures (magnetic field inversion or Feshbach resonances). This opens the possibility to study the physics of NT in a flexible experimental environment. For instance, the thermalization of two beams at different laser wavelengths interacting through the fiber nonlinearity can be exploited to achieve an efficient optical refrigeration: A highly incoherent speckled beam at NT can be cooled through its thermalization with a coherent beam towards a highly coherent state without any power loss. In contrast with usual beam cleaning at positive temperature where the energy is conserved, here the cooling process is featured by an energy transfer from the incoherent to the coherent beam, which significantly improves the gain of coherence of the incoherent beam.
Following this idea, one can explore the meaning of a *thermostat at NT* [10]: If the NT incoherent beam has a power much larger than the partially coherent beam, it will play the role of a NT thermal reservoir for such a partially coherent beam.

The versatile optical experimental environment proposed in this work also opens the possibility to study controversies about NTs, such as thermodynamic engines featured by Carnot cycles operating between temperatures of opposite signs, in relation with the generalized Kelvin-Planck formulation of the second law of thermodynamics stating that it is not possible to completely transform work into heat at NT [10].

*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).

# SUPPLEMENTARY MATERIAL

## EXPERIMENTAL SET-UP

The source is a Nd:YAG laser delivering subnanosecond pulses (400ps) at $\lambda = 1064$nm. The laser beam is passed through a spiral phase plate (Thorlabs) to generate a doughnut-like ring-shaped beam, and subsequently through a diffuser before injection of the speckle beam into the MMF, see Fig. 6. The diffuser plate is placed in the vicinity of the Fourier plane of a 4f-optical system. The near-field (NF) intensity distribution of the fiber output beam was magnified and imaged on a first CCD camera owing to a two lens telescope optical system, with $f_2 = 8$ mm and $f_3 = 150$ mm. 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 the magnified image was obtained by placing it in the object focal-plan of a lens $f_4 = 250$ mm and using a second CCD camera positioned in its image (Fourier) focal-plan.

We have computed analytically the propagation of the optical wave throughout the setup of our detection scheme, according to Fig. 6 (lower part). If $\psi_0(\mathbf{r})$ is the optical field amplitude at the fiber output ($\mathbf{r} = (x, y)$), then we have in the NF plane:

$$ \psi_{\text{NF}}(\mathbf{r}) = -\rho^{-1}\psi_0(-\mathbf{r}/\rho), $$

with $\rho = f_3/f_2$ the magnification factor. In the FF plane, the wave amplitude reads

$$ \psi_{\text{FF}}(\mathbf{u}) = \frac{i\rho}{\lambda f_4} \int d\mathbf{r} \psi_0(\mathbf{r}) \exp[-i2\pi(-\rho)\mathbf{r} \cdot \mathbf{u}/(\lambda f_4)], $$

<figure><img src="image_7.png" /><figcaption>FIG. 6: Setup. laser, optical isolator, half-wave plate and polarizer, lenses for magnification and imaging ($f_j$), spiral phase plate (V), diffuser (D), graded-index MMF, and cameras for near- and far-field detections (Cam).</figcaption></figure>

which corresponds to the Fourier transform of the field amplitude at the fiber output (note that the constant phase prefactor plays no role because the camera records the intensity). We note that: (i) The optical amplitude in the NF plane is an exact magnification of the wave amplitude at the output of the MMF; (ii) the optical amplitude in the FF detection plane exactly corresponds to the Fourier transform of the amplitude at the fiber output. Then, the experimental setup for the detection of the NF and FF intensities does not introduce detrimental spurious transverse phases profiles in the optical field, e.g., related to optical free propagation in air or phase shifts due to the presence of additional lenses.

## Multimode fiber

The refractive index profile of the graded-index MMF exhibits a parabolic shape in the fiber core with a maximum core index (at the center) of $n_{co} \approx 1.472$ and $n_{cl} \approx 1.457$ for the cladding at the pump wavelength of $1064$nm (fiber radius $R = 15\mu m$). The fiber length is $L = 12$ m. The trapping parabolic potential reads $V(\mathbf{r}) = q|\mathbf{r}|^2$ for $|\mathbf{r}| \le R$ and $q = k_0(n_{co}^2 - n_{cl}^2)/(2n_{co}R^2)$, $k_0 = 2\pi/\lambda$ the laser wave-number. The fundamental mode energy level is $\beta_0 = 2\sqrt{\alpha q}$, with $\alpha = 1/(2n_{co}k_0)$. The MMF guides $M = 45$ modes (i.e. $g = 9$ groups of degenerate modes).

The truncation of the potential introduces a frequency cut-off in the FF spectrum $k_c = (2\pi/\lambda)\sqrt{n_{co}^2 - n_{cl}^2}$ [23]. The conservation of $N$ and $E$ through propagation in the MMF (see Fig. 7) shows that the coupling from guided modes to leaky radiation modes in the cladding is negligible. This is not surprising as the efficiency of such a coupling can be shown to be very small. Indeed, the MMF has a core (radius $15\mu m$), a cladding (radius $62.5\mu m$), and a highly absorbing polymer-coating with refractive index larger than the core. Then leaky radiation modes in the cladding are rapidly absorbed during propagation due to their large penetration in the polymer-coating: we measured a typical absorption length $L_{abs}$ of $\approx 15$cm. Let us consider the coupled amplitude equations describing the four-wave mixing between four modes including a leaky mode. The four mode amplitudes $a_1, a_2, a_3, a_4$
(where $a_4$ stands for the leaky mode amplitude) satisfy Eqs.(10.2) in [44], in which we need to add an absorption term of the form $-a_4/L_{abs}$ in the equation (10.2.4) in [44] for $a_4$. Since absorption is strong the overdamped limit is valid and we get $a_4 = in_2k_0L_{abs}f_{4312}a_1a_2a_3^* \exp(-i\Delta\beta z)$ where $\Delta\beta = \beta_3 + \beta_4 - \beta_1 - \beta_2$, the parameter $n_2$ is the nonlinear-index coefficient, and $f_{4312}$ is an overlap integral between the four mode profiles. By substitution into one of the first three equations, say (10.2.1) in [44], we find that the mode amplitude $a_1$ experiences an effective absorption due to the coupling with the leaky mode that is given by $-4n_2^2k_0^2L_{abs}|f_{1234}|^2|a_2|^2|a_3|^2a_1$. This shows that this absorption is of the order of $-L_{abs}/L_{nl}^2a_1$ times a coefficient that is of the order of the square overlap integral between guided and leaky mode profiles. As the supports of these modes are very different (the guided modes are essentially supported in the core while the leaky modes are essentially supported in the cladding that is much larger), the square overlap integrals are small (smaller than the respective core-cladding ratio $(15/62.5)^4 \simeq 3 \times 10^{-3}$) and the effect of the coupling to the leaky modes onto the guided mode amplitudes can be neglected when $L_{nl} \simeq 20$cm and the propagation distance is $L = 12$m.

## Measurements of the energy $E$

From the measurements of the NF and FF intensity distributions, we have retrieved an accurate measurement of the power $N$ and the energy $E$ of the speckle beam. The NF intensity distribution $I_{NF}(\mathbf{r}) = |\psi|^2(\mathbf{r})$ provides a measurement of the power $N = \int I_{NF}(\mathbf{r})d\mathbf{r}$ and of the potential energy $E_{pot} = \int V(\mathbf{r})|\psi|^2(\mathbf{r})d\mathbf{r}$. The kinetic energy $E_{kin} = \alpha \int |\nabla\psi|^2(\mathbf{r})d\mathbf{r}$ is retrieved from the FF intensity distribution $I_{FF}(\mathbf{k}) = |\tilde{\psi}|^2(\mathbf{k})$. This provides the measurement of the (linear) energy (Hamiltonian) $E = E_{pot} + E_{kin}$.

## Conservation of $N$ and $E$ through propagation in the MMF for $E > E_*$ (negative temperature region)

Power conservation has been verified by keeping fixed the conditions of injection of the speckle beam into the MMF: We measured $N_{out}$ at $L = 12$m, and then $N_{in}$ by cutting the fiber at 20cm, and we always obtained $(N_{in} - N_{out})/N_{moy} < 1\%$. The experimental verification of energy conservation requires both the NF and FF intensity measurements. The NF and FF intensities are recorded at the fiber output at $L = 12$m, which gives $E_{out}$. Without altering the fiber launch conditions, the fiber is cut to 20cm to record the input NF and FF intensities, which gives $E_{in}$. The measurements of $E_{in}$ and $E_{out}$ then refer to an individual realization of the speckle beam (without average over the realizations). Fig. 7 shows that the conservation of the energy is well verified for $E > E_*$, i.e. in the negative temperature region. The

<figure><img src="image_8.png" /><figcaption>FIG. 7: Conservation of the energy through propagation in the MMF in the negative temperature region.</figcaption></figure>

(a) Measurements of the energy at the input of the MMF (blue triangles), and at the output of the MMF (red triangles): The energy $E/N$ is conserved through the propagation in the MMF over a broad range of variation of $E/E_{min}$. We recall that negative temperatures $T < 0$ occur for $E > E_*$ with $E_*/E_{min} \simeq 6.33$, see Fig. 1.

energy $E$ is varied owing to the diffuser before injection into the MMF, see Fig. 6.

## MODAL DECOMPOSITION

### Phase retrieval

The procedure of mode decomposition is based on the well-known Gerchberg-Saxton algorithm [41, 50, 51]. From the measurements of the NF and FF intensity distributions in the experiment, it allows us to retrieve the transverse phase profile of the field. The resulting complex field is subsequently projected onto the fiber modes, to get the complete modal distribution of the experimental optical beam. The algorithm is known to be accurate although it is not efficient in terms of computational cost. Indeed it is a local search algorithm that updates iteratively the unknown phase profile of the field and it is usually necessary to consider several initial phase guesses. We have, therefore, carried out a detailed preliminary analysis of the algorithm by performing numerical simulations that reproduce our experimental configuration in order to prove that the phase retrieval and modal distribution estimation are reliable.

### Error introduced by the Gerchberg-Saxton algorithm

To evaluate accurately the error in the phase-retrieval algorithm, we have reproduced in detail the experimental procedure as follows:

i) We consider a particular value of the energy $E$ ($E > E_*$ so that $T < 0$). Throughout the procedure the power $N$ is set constant. The pair $(E, N)$ determines uniquely $(T, \mu)$ and thus the exact RJ distribution at equilibrium $n_p^{RJ} = T/(\beta_p - \mu)$.
<figure><img src="image_9.png" /><figcaption>FIG. 8: **Phase retrieval.** Example of numerically generated near-field intensity distribution (a), and its reconstruction (b). Original phase field (c) and the corresponding phase field reconstructed from the Gerchberg-Saxton algorithm (d).</figcaption></figure>

ii) We generate from $n_p^{\text{RJ}}$ a realization of speckle beam $\psi(\mathbf{r}) = \sum_p a_p u_p(\mathbf{r})$, where $a_p$ is a complex Gaussian random variable with variance $\langle |a_p|^2 \rangle = n_p^{\text{RJ}}$ ($a_p = a_p^{(r)} + ia_p^{(i)}$ with $a_p^{(r)}$ and $a_p^{(i)}$ real Gaussian independent random variables with mean zero and variance $n_p^{\text{RJ}}/2$). We recall that $u_p(\mathbf{r})$ are the fiber modes. Then $\psi(\mathbf{r})$ is a particular realization of a complex speckle field at exact RJ equilibrium.

iii) The particular realization $\psi(\mathbf{r})$ is highly resolved numerically. We mimic the impact of the finite resolution of the camera used in the experiment. From $\psi(\mathbf{r})$ we compute the NF and FF intensity distributions $|\psi(\mathbf{r})|^2$ and $|\hat{\psi}(\mathbf{k})|^2$. We sample the NF and FF intensity distributions with the finite number of points available in the camera (1024²) in r-space and k-space and ≈ 950 points for the dynamics range in intensity. We apply the Gerchberg-Saxton algorithm to retrieve the sampled phase profile. Due to the errors introduced by the sampling of the camera and by the phase-retrieval algorithm, the resulting complex field $\psi^{\text{exp}}(\mathbf{r})$ may differ from the generated speckle beam $\psi(\mathbf{r})$.

We report in Fig. 8 the numerically generated near-field intensity distribution in one numerical simulation (a) and its reconstruction (b), the original phase profile (c) and the reconstructed phase profile (d). It is clear that the reconstruction is very good. We will see below more quantitatively that the error is indeed negligible.

iv) We project the complex field $\psi^{\text{exp}}(\mathbf{r})$ onto the fiber modes to get the complex modal coefficient $a_p^{\text{exp}}$ and distribution $|a_p^{\text{exp}}|^2$. Due to the errors introduced by the sampling of the camera and by the phase-retrieval algorithm, this modal distribution may differ from the mode

<figure><img src="image_10.png" /><figcaption>FIG. 9: **Error in the modal decomposition.** The mode decomposition is based on the Gerchberg-Saxton algorithm, whose error is quantified by the distance $D_{\text{err}}^Q$ to the exact distribution, see Eq.(4). (a) $D_{\text{err}}^Q$ vs the energy $E$ (for $Q = 50$ realizations), by accounting for the sampling due to the limited resolution of the camera (red), and without the sampling (blue). Note in (a) that an increase of the randomness of the speckle beam (i.e., increase of $E$) does not increase the error. (b) $D_{\text{err}}^Q$ vs number of realizations $Q$ (for $E/E_{\text{min}} = 7$): The error decreases with the number $Q$ of realizations of the speckles. The green line reports the theoretical estimate of the error given in Eq.(6). Note that $D_{\text{err}}^Q$ is bounded, $0 \le D_{\text{err}}^Q \le 1$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10: **Near-field and far-field experimental intensities, and corresponding reconstructed intensity distributions.** Near-field (NF) intensity recorded in the experiment for a single realization of a speckle (a), and corresponding far-field (FF) intensity distribution (c). Corresponding NF intensity (b), and FF intensity (d), reconstructed from the Gerchberg-Saxton algorithm.</figcaption></figure>
distribution $|a_p|^2$ used to generate the speckle beam $\psi(\mathbf{r})$. As we will show below, this error is negligible.

v) We repeat the steps ii)-iv) Q times, each with a different realization of the speckle beam (i.e., with a different realization $a_p^j$ of $a_p$, for $j = 1, \dots, Q$). The procedure then gives Q distributions $|a_p^{\exp,j}|^2$, $j = 1, \dots, Q$. We compute the empirical averages $n_p^{\exp,Q} = (1/Q) \sum_{j=1}^Q |a_p^{\exp,j}|^2$. We anticipate that, for Q large enough, these empirical averages should be close to the theoretical values $n_p^{\text{RJ}}$. We introduce the estimation error:

$$ \mathcal{D}_{\text{err}}^{Q} = \frac{\sum_{p} |n_{p}^{\text{exp},Q} - n_{p}^{\text{RJ}}|}{\sum_{p} n_{p}^{\text{exp},Q} + n_{p}^{\text{RJ}}}. \quad (4) $$

Let us imagine for a while that the phase-retrieval algorithm is perfect and that the sampling error due to the camera is absent. Then, for each realization $j = 1, \dots, Q$, we have $|a_p^{\exp,j}|^2 = |a_p^j|^2$ exactly. Thus, the random variables $|a_p^{\exp,j}|^2$ are independent and follow exponential distributions with mean $n_p^{\text{RJ}}$. Consequently, the empirical quantities $Z_p^Q = n_p^{\exp,Q}/n_p^{\text{RJ}}$ are independent and identically distributed with the gamma probability distribution $\Gamma(Q, Q)$ (the law of the sum of Q independent variables with exponential distribution and mean $1/Q$) and the estimation error is

$$ \mathcal{D}_{\text{err}}^{Q} = \frac{\sum_{p} n_{p}^{\text{RJ}} |Z_{p}^{Q} - 1|}{\sum_{p} n_{p}^{\text{RJ}} (Z_{p}^{Q} + 1)}, \quad (5) $$

which gives when the number of modes is large enough ($Z^Q$ follows the $\Gamma(Q, Q)$ distribution):

$$ \mathcal{D}_{\text{err}}^{Q} \approx \frac{\mathbb{E}[|Z^Q - 1|]}{\mathbb{E}[Z^Q + 1]} = \frac{Q^{Q-1}}{(Q-1)!} e^{-Q}. \quad (6) $$

For $Q \ge 8$, we have $\mathcal{D}_{\text{err}}^Q \approx 1/\sqrt{2\pi Q}$.

We have carried out numerical simulations with our implementation of the phase-retrieval algorithm (using multiple initial phase guesses) and with the sampling error of the camera. The results of the distance $\mathcal{D}_{\text{err}}^Q$ vs energy $E$ are reported in Fig. 9 with different numbers $Q$ of realizations per energy. We can see in panel (b) of Fig. 9 that the errors correspond to the theoretical error Eq.(6) when the phase-retrieval algorithm makes no error.

The error introduced by the Gerchberg-Saxton algorithm has been computed by increasing the amount of complexity in the speckle pattern, i.e., by increasing the energy $E$. We can see in Fig. 9(a) that the error does not increase when the energy $E$ increases.

In the experiments we have typically 35 to 70 independent realizations of speckle beams for a given small energy interval $[E - \Delta E, E + \Delta E]$ with $\Delta E = 0.125 E_{\text{min}}$, so we can expect that the errors (due to the phase retrieval algorithm and the camera sampling) are small. Error bars with relative standard deviations of the order of $1/\sqrt{2\pi Q} \sim 6\%$ could be added in Fig. 2 but they are too small to be visible.

<figure><img src="image_12.png" /><figcaption>FIG. 11: Experimental attraction to NT RJ equilibrium. Distance $\mathcal{D}_{\text{RJ}}$ [defined in Eq.(7)] to the RJ equilibrium distribution computed from the experimental data at the fiber input (blue), and fiber output (red), for different values of the energy $E$. The significant reduction of $\mathcal{D}_{\text{RJ}}$ from input to output measurements shows the attraction to the RJ equilibrium for $T < 0$. Note that $\mathcal{D}_{\text{RJ}}$ in Eq.(7) is bounded, $0 \le \mathcal{D}_{\text{RJ}} \le 1$. We recall that negative temperatures $T < 0$ occur for $E > E_*$ with $E_*/E_{\text{min}} \approx 6.33$, see Fig. 1.</figcaption></figure>

To complete our study, we report in Fig. 10 the near-field and far-field intensity distributions recorded during one of the experiments (left plots) and the corresponding reconstructed intensities from the Gerchberg-Saxton algorithm (right plots).

## Experimental convergence to the NT RJ distribution

We have quantified in our experimental results the attraction to the NT equilibrium by using Eq.(4), which provides a ‘distance’ to the RJ distribution

$$ \mathcal{D}_{\text{RJ}} = \frac{\sum_{p} |n_{p}^{\text{exp}} - n_{p}^{\text{RJ}}|}{\sum_{p} n_{p}^{\text{exp}} + n_{p}^{\text{RJ}}}. \quad (7) $$

We report in Fig. 11 the distance $\mathcal{D}_{\text{RJ}}$ computed for the experimental data averaged over the realizations $n_p^{\text{exp}}$ at the fiber input (blue), and the fiber output (red), for different energies $E$. The strong reduction of the distance $\mathcal{D}_{\text{RJ}}$ from the fiber input to the output confirms the process of NT thermalization, which is demonstrated over a broad range of values of the energy $E$.

## Polarization effects

The polarization state of the optical beam changes as it propagates through the MMF. The field at the output of the MMF is projected onto a basis of orthogonal linear polarizations. The corresponding NF and FF intensity distributions are recorded along the orthogonal linear polarizations. For each polarization, we apply the mode decomposition procedure based on the Gerchberg-Saxton presented above. In this way, we retrieve the transverse phase profile of the field for the two orthogonal polarizations. The complex field along each polarization is subsequently projected over the fiber modes, to get the modal populations for each polarization, $n_p^{(x)}$ and $n_p^{(y)}$. The modal distribution is obtained
<figure><img src="image_13.png" /><figcaption>FIG. 12: Polarization effects in modal decomposition. Experimental modal distribution retrieved from the Gerchberg-Saxton algorithm, without splitting the orthogonal polarization (blue line), by splitting the orthogonal polarization at the fiber output (dashed red line). Note that a single realization of the speckle beam has been considered, which explains the discrepancy with the RJ equilibrium distribution (dashed green line).</figcaption></figure>

by summing the contributions of the two polarizations, $n_p^{\text{pol}} = (N_x n_p^{(x)} + N_y n_p^{(y)})/(N_x + N_y)$, where $N_{x,y}$ denote the power along the two polarizations. We compare in Fig. 12 the modal distribution retrieved by following this procedure ($n_p^{\text{pol}}$), with the modal distribution ($n_p$) retrieved without separating the polarization states of the field. This comparison is reported in Fig. 12 for the same (single) realization of the speckle beam. We can remark in Fig. 12 that the two distributions are almost identical. In all cases analyzed, we have always observed the same good agreement. Given the large number of realizations of speckle beams ($\simeq 300$) recorded and analyzed in our experiments, we didn't perform the polarization modal decomposition.

## ANALOGUE CONDENSATION AT NT

### Thermodynamics relations

We define the thermodynamic relations used to plot Fig. 1. We consider a field at equilibrium with the RJ distribution $n_p^{\text{RJ}} = T/(\beta_p - \mu)$, with $N = \sum_p n_p^{\text{RJ}}$ and $E = \sum_p \beta_p n_p^{\text{RJ}}$, where $\beta_p = \beta_{p_x, p_y} = \beta_0(p_x + p_y + 1)$ are the eigenvalues of the truncated parabolic potential ($\beta_p \le \beta_{\text{max}}$), and the index $\{p\}$ labels the two integers $(p_x, p_y)$ that specify a mode. Here and below, the sum over the modes $\sum_p$ is carried over the set $0 \le p_x + p_y < g$, where $g = \beta_{\text{max}}/\beta_0$ is the number of groups of non-degenerate modes, with $M = g(g+1)/2$ the total number of modes. Because of the constraint $n_p^{\text{RJ}} = T/(\beta_p - \mu) > 0$, a NT equilibrium state $T < 0$ requires that $\mu > \beta_{\text{max}}$.

We start from the equilibrium entropy $\tilde{S}^{\text{eq}} = \sum_p \log(n_p^{\text{eq}})$ – note that at equilibrium it coincides with the nonequilibrium entropy verifying the H theorem of entropy growth. It proves convenient to shift the entropy by a constant $S^{\text{eq}} = \tilde{S}^{\text{eq}} - M \log N$, so that by

<figure><img src="image_14.png" /><figcaption>FIG. 13: Analogue effect of condensation at NTs: (a) Chemical potential $\mu/\beta_0 - 1$ vs energy $E/E_{\text{min}}$ for the MMF used in the experiments ($g = 9$), from Eq. (9). Note the asymptotic behaviors $\mu \to \beta_0^-$ for $E \to E_{\text{min}}$, and $\mu \to \beta_{\text{max}}^+$ for $E \to E_{\text{max}}$, which lead respectively to the macroscopic population of the lowest energy level ($\beta_0$), and highest energy level ($\beta_{\text{max}}$). The horizontal dashed line denotes $\mu = \beta_{\text{max}}$ and the vertical one $E = E_{\text{*}}$. (b) Condensate fraction in the highest energy level $\tilde{n}_g^{\text{RJ}}/N$ vs energy $E$: As the energy increases $E \to E_{\text{max}}$ (or equivalently the temperature increases $T \to 0^-$), the condensate fraction increases $\tilde{n}_g^{\text{RJ}}/N \to 1$. The curves are obtained from Eqs. (14-15), see the text. The critical behavior of the transition to condensation becomes apparent by increasing the number of modes $g$.</figcaption></figure>

using $T = N/\sum_p (\beta_p - \mu)^{-1}$, we can write

$$S(\mu) = -\sum_p \log(\beta_p - \mu) - M \log \left( \sum_p \frac{1}{\beta_p - \mu} \right) \quad (8)$$

$$\frac{E(\mu)}{E_{\text{min}}} = \frac{\sum_p \frac{\beta_p}{\beta_p - \mu}}{\sum_p \frac{\beta_0}{\beta_p - \mu}} \quad (9)$$

$$\frac{T(\mu)}{E_{\text{min}}} = \frac{1}{\sum_p \frac{\beta_0}{\beta_p - \mu}} \quad (10)$$

The parametric plot with resp. to $\mu$ of (8) and (9) gives $S(E)$ in Fig. 1(b); the corresponding parametric plot of (9) and (10) gives $T$ vs $E$ in Fig. 1(c).

## NT states in the thermodynamic limit

The thermodynamic limit is defined by $N \to \infty$, $\beta_0 \to 0$ with $N\beta_0^2 = \text{const}$ and $\beta_{\text{max}} = \text{const}$. In this limit the discrete sums over the modes are replaced by continuous integrals, namely $N = \sum_p n_p^{\text{eq}} \to$
$(T/\beta_0^2) \int_0^{\beta_{\max}} dx \int_0^{\beta_{\max}-x} dy (x+y+\beta_0-\mu)^{-1}$, which gives

$$N\beta_0^2 = T\beta_{\max}\left(1 - z \log\left(\frac{z}{z-1}\right)\right), \quad (11)$$

where $z = \mu/\beta_{\max}$. A negative temperature equilibrium state ($T < 0$) is characterized by $z > 1$. It can exist in the thermodynamic limit with the constraint $N\beta_0^2 = \text{const} > 0$ provided that $1 - z \log (z/(z-1)) < 0$. This inequality is always verified for $z > 1$, which confirms the existence of NT equilibrium states in the thermodynamic limit.

## NT condensation in the highest energy level

Positive temperature $T > 0$. We start by briefly summarizing the usual positive temperature condensation in the lowest energy level. This effect of condensation originates in the singularity of the RJ distribution. Indeed, the denominator of $n_p^{\text{RJ}} = T/(\beta_p - \mu)$ vanishes for the lowest energy level when $\mu = \beta_0$. This gives rise to an analogue effect of condensation: As the energy decreases below a critical value $E_{\text{crit}}$ (or $T < T_{\text{crit}}$), then $\mu \to \beta_0^-$ (see Fig. 13(a)) and the singular behavior of the RJ distribution is regularized by the macroscopic population of the fundamental mode:

$$n_0^{\text{RJ}}/N \to 1 \quad \text{as} \quad E \to E_{\text{min}} \quad (\text{or } T \to 0^+). \quad (12)$$

It has been shown that this condensation-like effect is a phase transition that occurs in the thermodynamic limit, see Ref. [23].

Negative temperature, $T < 0$. In the NT region, we reveal an inverted condensation-like effect, which is characterized by a macroscopic population of the highest energy level. Aside from the singularity for $\mu = \beta_0$ discussed here above for $T > 0$, the RJ distribution $n_p^{\text{RJ}} = T/(\beta_p - \mu)$ also exhibits a singularity (vanishing denominator) for the highest energy level when $\mu = \beta_{\max}$ (see Fig. 13(a)). We have shown in Fig. 5, that the highest energy level $g = 9$ becomes macroscopically populated as the energy increases:

$$\tilde{n}_g/N \to 1 \quad \text{as} \quad E \to E_{\text{max}} \quad (\text{or } T \to 0^-), \quad (13)$$

with $\tilde{n}_g = \sum_{p,p_x+p_y+1=9} n_p$. This condensation-like effect does not occur in the thermodynamic limit. We pose $\mu = \beta_{\max} + \varepsilon$, with $\varepsilon > 0$. Eq. (11) can be written in the limit $\varepsilon \to 0^+$: $N\beta_0^2 \simeq T\beta_{\max}(1 - \log(\beta_{\max}/\varepsilon))$. By keeping $N\beta_0^2 = \text{const}$, the chemical potential $\mu$ reaches $\beta_{\max}^+$ for a vanishing temperature $T \to 0^-$, i.e., condensation does not occur in the thermodynamic limit. However, an analogue effect of condensation occurs through a macroscopic population of the highest energy level. By setting $\beta_p = \beta_{\max}$ in the RJ distribution, then $n_p^{\text{RJ}} = -T/(\mu - \beta_{\max})$ denotes the power in one mode of the highest energy level. The total power in the highest energy level ($g$-fold degenerate) $\tilde{n}_g^{\text{RJ}} = gn_p^{\text{RJ}}$ then

<figure><img src="image_15.png" /><figcaption>FIG. 14: Simulations of NLS equation and wave turbulence kinetic equation: Simulations of the NLS equation (dashed lines), and wave turbulence kinetic equation (solid lines) showing the evolutions during the propagation (in $z$) of the power within each of the $g=9$ groups of degenerate modes of the fiber that has been used in the experiments, for $E/E_{\text{min}} = 8$ (a), $E/E_{\text{min}} = 7$ (b). The horizontal dashed black line denotes the population for the higher mode group at complete RJ thermal equilibrium. See the text for details on initial conditions and parameters.</figcaption></figure>

reads

$$\frac{\tilde{n}_g^{\text{RJ}}}{N}(\mu) = \frac{g}{(\mu - \beta_{\max}) \sum_p (\mu - \beta_p)^{-1}}, \quad (14)$$

$$\frac{E(\mu)}{E_{\text{min}}} = \frac{\sum_p \frac{\beta_p}{\beta_p - \mu}}{\sum_p \frac{\beta_0}{\beta_p - \mu}}. \quad (15)$$

The parametric plot of (14) and (15) with respect to $\mu$ provides the condensate fraction reported in Fig. 5. Fig. 13(b) shows the condensate fraction, $\tilde{n}_g^{\text{RJ}}/N$ vs $E$, by increasing the number of modes, i.e., by approaching the thermodynamic limit. The critical behavior of the condensation curve looks similar to that of a phase transition, thought strictly speaking phase transitions only occur in the thermodynamic limit. Nevertheless, if one takes the macroscopic occupation of an energy level as the essential characteristic of condensation, then NTs are characterized by a condensation-like effect into the highest-energy level.

## WAVE TURBULENCE SIMULATIONS

The starting point is the NLS equation governing light propagation in MMFs [39, 42]. By expanding the random wave into the fiber modes and considering the dominant contribution of weak disorder, the polarized modal components $a_p = (a_{p,x}, a_{p,y})^T$ are governed by [39]:

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

where $\mathbf{P}_p(a) = \sum_{l,m,n} S_{plmn} \left( \frac{1}{3} a_l^T a_m a_n^* + \frac{2}{3} a_n^\dagger a_m a_l \right)$, $S_{plmn}$ denoting the spatial overlap among the modes. The introduction of disorder is important in order to avoid strong phase correlations among the modes that are related to Fermi-Pasta-Ulam recurrences [43], which inhibit the thermalization process [38]. We consider the
most general form of disorder that conserves the power: The Hermitian matrices $D_p(z)$ are expanded into the Pauli matrices that 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$, where $\sigma_j$ ($j = 1, 2, 3$) are the Pauli matrices ($\sigma_0$ the identity matrix), while $\nu_{p,j}(z)$ are independent and identically distributed real-valued random processes, with variance $\sigma^2$ and correlation length $\ell_c$. The corresponding characteristic length scale of disorder is $L_d = 1/\Delta\beta$, with $\Delta\beta = \sigma^2\ell_c$, see Ref.[39] for details.

The nonequilibrium process of NT thermalization can be described by a wave-turbulence kinetic equation. The kinetic equation was derived from the modal NLS Eq.(16) in the weakly nonlinear regime of the experiment, $L_{lin} \sim 1/\beta_0 \ll L_{nl} \sim 1/(\gamma N)$. It describes the nonequilibrium evolution of the averaged modal components $n_p(z) = \langle |\mathbf{a}_p|^2(z) \rangle$ [39]:

$$ \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 (17) $$

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$ and $\Delta\omega_{lp} = \beta_l - \beta_p$. The term $\delta^K(\Delta\omega_{lmnp})$ denotes the four-wave frequency resonance $\Delta\omega_{lmnp} = \beta_l + \beta_m - \beta_n - \beta_p$, with $\delta^K(\Delta\omega_{lmnp}) = 1$ if $\Delta\omega_{lmnp} = 0$, and zero otherwise. The kinetic Eq.(17) conserves $N = \sum_p n_p(z)$, the energy $E = \sum_p \beta_p n_p(z)$ and exhibits a $H$-theorem of entropy growth, $\partial_z S_{kin}(z) \ge 0$, for the nonequilibrium entropy $S_{kin}(z) = \sum_p \log(n_p(z))$. Accordingly, it describes an irreversible evolution of the speckle beam to the RJ equilibrium distribution realizing the maximum of entropy, $n_p^{RJ} = T/(\beta_p - \mu)$ [39].

We have considered in the simulations the MMF used in the experiments, see Sec. . The initial condition consists of a speckle beam whose correlation length is varied in such a way to fix a desired value of the energy $E$. The considered parameters are $\ell_c = 0.3$ m and $2\pi/\sigma = 2.1$ m [39]. The results of the simulations of the NLS equation and the corresponding simulations of the wave turbulence kinetic equation are reported in Fig. 14. They show the process of thermalization to the negative temperature RJ equilibrium state predicted by the theory. The simulations qualitatively reproduce the experimental results, although a power of 20kW has been considered to accelerate the dynamics. The purely spatial model considered here then captures the essential features of the condensation process reported experimentally.

At variance with NLS simulations that are stochastic and thus exhibit fluctuations, the simulations of the wave turbulence kinetic equation are deterministic (free of fluctuations). A good agreement between NLS and kinetic simulations is obtained without using any adjustable parameter. The wave turbulence kinetic equation then provides a nonequilibrium description of the

process of negative temperature thermalization observed experimentally.

[1] L. Onsager, Statistical hydrodynamics, Il Nuovo Cimento (1943-1954) 6, 279 (1949).

[2] N. F. Ramsey, Thermodynamics and statistical mechanics at negative absolute temperatures, Physical Review 103, 20 (1956).

[3] J. Dunkel and S. Hilbert, Consistent thermostatistics forbids negative absolute temperatures, Nature Physics 10, 67 (2014).

[4] S. Calabrese and A. Porporato, Origin of negative temperatures in systems interacting with external fields, Physics Letters A 383, 2153 (2019).

[5] D. Frenkel and P. B. Warren, Gibbs, Boltzmann, and negative temperatures, American Journal of Physics 83, 163 (2015).

[6] P. Buonsante, R. Franzosi, and A. Smerzi, On the dispute between Boltzmann and Gibbs entropy, Annals of Physics 375, 414 (2016).

[7] A. Puglisi, A. Sarracino, and A. Vulpiani, Temperature in and out of equilibrium: A review of concepts, tools and attempts, Physics Reports 709, 1 (2017).

[8] L. Cerino, A. Puglisi, and A. Vulpiani, A consistent description of fluctuations requires negative temperatures, Journal of Statistical Mechanics: Theory and Experiment 2015, P12002 (2015).

[9] E. Abraham and O. Penrose, Physics of negative absolute temperatures, Physical Review E 95, 012125 (2017).

[10] M. Baldovin, S. Iubini, R. Livi, and A. Vulpiani, Statistical mechanics of systems with negative temperature, Physics Reports 923, 1 (2021).

[11] M. Onorato, G. Dematteis, D. Proment, A. Pezzi, M. Ballarin, and L. Rondoni, Equilibrium and nonequilibrium description of negative temperature states in a one-dimensional lattice using a wave kinetic approach, Phys. Rev. E 105, 014206 (2022).

[12] Note that, particles in a laser have an inverted energy population which, however, only exists as long as the laser is continuously pumped (by switching off the pump, all atoms decay back into the lower state and their energy goes into the laser beam). The laser medium is in a steady state, but not in thermal equilibrium and therefore cannot have a temperature.

[13] E. M. Purcell and R. V. Pound, A nuclear spin system at negative temperature, Physical Review 81, 279 (1951).

[14] S. Braun, J. P. Ronzheimer, M. Schreiber, S. S. Hodgman, T. Rom, I. Bloch, and U. Schneider, Negative absolute temperature for motional degrees of freedom, Science 339, 52 (2013).

[15] S.P. Johnstone, A.J. Groszek, P.T. Starkey, C.J. Billington, T.P. Simula, K. Helmerson, Evolution of large-scale flow from turbulence in a two-dimensional superfluid, Science 364, 1267 (2019).

[16] G. Gauthier, M. T. Reeves, X. Yu, A. S. Bradley, M.A. Baker, T.A. Bell, H. Rubinsztein-Dunlop, M.J. Davis, T.W. Neely, Giant vortex clusters in a two-dimensional quantum fluid, Science 364, 1264 (2019).

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

[18] 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).

[19] C. Conti, M. Leonetti, A. Fratalocchi, L. Angelani, G. Ruocco, Condensation in disordered lasers: Theory, 3d+1 simulations, and experiments, Phys. Rev. Lett. **101**, 143901 (2008).

[20] J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz, Bose-Einstein condensation of photons in an optical microcavity, Nature (London) **468**, 545 (2010).

[21] R. Weill, A. Bekker, B. Levit, and B. Fischer, Bose-Einstein condensation of photons in an erbium-ytterbium co-doped fiber cavity, Nature Communications **10**, 747 (2019).

[22] J. Bloch, I. Carusotto, M. Wouters, Spontaneous coherence in spatially extended photonic systems: Non-Equilibrium Bose-Einstein condensation, Nature Physics Reviews **4**, 470 (2022).

[23] 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).

[24] H. Pourbeyram, P. Sidorenko, F. Wu, L. Wright, D. Christodoulides, F. Wise, Direct measurement of thermalization to Rayleigh-Jeans distribution in optical beam self cleaning, Nature Physics **18**, 685 (2022).

[25] 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).

[26] 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).

[27] 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).

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

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

[30] 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).

[31] 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).

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

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

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

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

[36] 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 (2014).

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

[38] 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).

[39] 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).

[40] 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).

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

[42] 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).

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

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

[45] B. Rumpf, Stable and metastable states and the formation and destruction of breathers in the discrete nonlinear Schrödinger equation, Physica D **238**, 2067 (2009).

[46] S. Iubini, L. Chirondojan, G.-L. Oppo, A. Politi, P. Politi, Dynamical Freezing of Relaxation to Equilibrium, Phys. Rev. Lett. **122**, 084102, (2019).

[47] N.K. Efremidis, D.N. Christodoulides, Fundamental entropic processes in the theory of optical thermodynamics, Phys. Rev. A **103**, 043517 (2021).

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

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

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

[51] 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).