2023

[1]Mario Zitelli

[1]Department of Information Engineering, Electronics and Telecommunications,
Università degli Studi di Roma Sapienza, Via Eudossiana 18, Rome, 00184, RM,
Italy

# Statistics of modal condensation in nonlinear multimode fibers

<EMAIL_ADDRESS>Fabio Mangini<EMAIL_ADDRESS>Stefan
Wabnitz<EMAIL_ADDRESS>*

###### Abstract

Optical pulses propagating in multimode optical fibers are affected by linear
disorder and nonlinearity, and experience chaotic exchange of power among
modes. On the other hand, complex systems can attain steady states
characterized by energy condensation into single as well multiple sub-systems.
In this work, we study beam propagation in multimode optical fibers in the
presence of linear random mode coupling and Kerr nonlinearity; both effects
lead to a mode power redistribution at the fiber output. We use a new 3D mode
decomposition method to obtain, with unprecedented accuracy, measurements of
the modal distribution from long spans of graded-index fiber; we perform
numerical simulations using a new model for the linear disorder; we introduce
a weighted Bose-Einstein law and show that it is suitable for describing
steady-state modal power distributions both in the linear and nonlinear
regimes. We show that, at power levels intermediate between the linear and the
soliton regimes, energy condensation is attained locally by the second, third
and fourth modal groups, before global condensation to the fundamental mode is
reached in the soliton regime. Our results extend the thermodynamic approach
to multimode fibers to unexplored optical states, which acquire the
characteristics of optical glass.

###### keywords:

multimode fibers, disorder, solitons, nonlinear optics

##  1 Introduction

Multimode (MM) optical fibers 6767859 ;
https://doi.org/10.1002/j.1538-7305.1972.tb01911.x ; Hasegawa:80  have
regained considerable interest over the last decade, motivated by the
potential for increasing the transmission capacity using the technique of
spatial-division multiplexing (SDM) Richardson-NatPhot-2013-94-2013 ;
Winzer:18 , and the possibility of up-scaling the pulse energy delivered from
fiber lasers wright2017spatiotemporal .

In recent years, a thermodynamic interpretation of beam propagation in
multimode optical fibers and systems has been formulated Wu2019-nn ;
Fusaro_PhysRevLett.122.123902 . In the specific example of a graded-index
(GRIN) MM fiber, photons are divided into indistinguishable energy packets:
statistical mechanics permits to estimate the equilibrium distribution of
these packets among degenerate groups of modes, with the same propagation
constant $\beta_{j}$ (m-1); such equilibrium distribution corresponds to
extrema of the entropy $S$ of the mode population.

On the other hand, it is well known that optical pulses propagating in a MM
fiber in the linear regime are affected by random-mode coupling (RMC)
Gloge:6774107 ; Marcuse1973LossesAI ; Ho:14 , which is induced by fiber
imperfections such as micro and macro-bending. Pulses in modes within
degenerate (i.e., with the same or nearly the same propagation constant)
groups separate in time from pulses in other groups, as a consequence of
inter-modal dispersion. 3D bullets carried by each group broaden in time
because of intra-modal or chromatic dispersion (Fig.1a). In the linear regime,
RMC effects can be modelled by power-flow equations Savovi2019PowerFI ;
Gloge:6774107 ; Olshansky:75 , which predict a diffusion of energy from
intermediate groups into both lower-order and high-order modes. These models
lead, after few hundreds meters of propagation, to a steady-state mode power
distribution showing a characteristic decrease of power as the mode order
grows larger.

At higher pulse energy, in MM fibers with anomalous chromatic dispersion,
optical solitons form with a specific pulsewidth $T_{\text{FWHM}}$ depending
on the pulse wavelength Zitelli2021a , and with energy given by
$E_{s}=1.76\lambda\lvert\beta_{2}(\lambda)\rvert
w_{e}^{2}/(n_{2}T_{\text{FWHM}})$, with $\beta_{2}(\lambda)$ (s2/m) the
chromatic dispersion, $n_{2}$ m2/W the Kerr nonlinear coefficient, and $w_{e}$
the modal effective waist. At telecom wavelength $\lambda=1550$ nm, the
soliton pulsewidth is $T_{\text{FWHM}}=120$ fs. Pulses corresponding to the
different modal groups get shorter and increase their peak power; Kerr and
Raman nonlinearities induce a slow but monotonic transfer of energy from each
group to the fundamental mode Zitelli2021 , leading to global condensation in
the ground state. When pulses separate in time, the mode coupling process is
eventually continued by the interplay of RMC and IM-FWM. After hundreds of
meters of propagation, a train of fundamental soliton bullets remains, which
experience Raman soliton self-frequency shift Gordon1986 ; Zitelli_9887813
(Fig.1c). The process can be viewed as a new fission mechanism, mediated by
the modal dispersion, as it will be illustrated in Sec 4.2.

At intermediate energy, 20% to 80% of the soliton $E_{s}$, a quasi-soliton
forms, characterized by pulses carrying the individual modal groups which are
initially overlapped in time, and reduce their pulsewidth to values
approaching the soliton $T_{\text{FWHM}}$. Mode power exchange is dominated by
the IM-FWM, while Raman self-frequency shift is not able to affect the pulse
wavelength yet; a conservative wave propagation system can still be assumed.
In this regime, local condensation of energy among the lower-order groups is
observed (Fig.1b); the modes assume a characteristic distribution, which
cannot be traced back to what has been observed in the linear or soliton
regimes. The local condensation observed in MM shows similarities to those in
disordered lasers PhysRevLett.101.143901  and Bose-Einstein condensation
Klaers_2011 .

In the following sections, a new weighted Bose-Einstein (BE) equation will be
demonstrated, which corresponds to extrema of the entropy of an optical
multimode system (Sec. 2). Experiments will be shown in long spans of graded-
index (GRIN) multimode fibers, from linear to soliton regimes (Sec. 3); modal
content emerging from the fiber will be analyzed using the weighted BE in
order to check the achievement of local or global condensation states.
Numerical power-flow simulations will be performed in the linear regime in
Sec. 4.1, and analyzed using the weighted BE. Finally, numerical simulations
using coupled-mode generalized nonlinear Schröedinger equations (GNLSE),
including nonlinearity and a new model for RMC (Sec. 4.2), will be used for
comparison with the experimental data. New steady states will be found to form
in the quasi-soliton regime, characterized by local condensates which can be
viewed as glassy states parisi_urbani_zamponi_2020 ; PhysRevB.83.134204 ; this
condition appears as intermediate between the low-energy disordered state, and
the high-energy, highly condensed soliton state.

Figure 1: Optical bullets forming after long spans of GRIN fiber in a) linear,
b) nonlinear and c) soliton regimes. A thermodynamic representation of the
three states is artistically represented.

##  2 Weighted Bose-Einstein Law

We start considering an optical multimode system including $Q$ groups of
degenerate modes, distributed over $g_{j}/2$ modes and 2 polarizations
($g_{j}$ is the group degeneracy), with $j=1,2,..,Q$; in the special case of a
graded-index multimode fiber (GRIN), it is $g_{j}=2,4,6,..,2Q$. Let $n_{j}$ be
the population of energy packets into the j-th group over 2 polarizations;
packets are not necessarily single photons ( groups of indistinguishable
photons can compose a single packet); hence, the value of $n_{j}$ is of the
order of 105 to 109 in practical applications. Extremization of the system’s
entropy $S$, with the proper choice of Lagrange multipliers, brings to the
weighted Bose-Einstein distribution for the modal power fractions (see Methods
Sec. 6.1)

$\lvert
f_{i}\rvert^{2}=\frac{2(g_{i}-1)}{g_{i}\gamma}\frac{1}{\exp\big{(}-\frac{\mu^{\prime}+\epsilon_{i}}{T}\big{)}-1}.$
(1)

In Eq. 1, $\lvert f_{i}\rvert^{2}=2n_{j}/(\gamma n_{0}g_{i})$ is the mean
modal power fraction over two polarizations, with $\gamma$ given by the total
number of energy packet $N=\gamma n_{0}$ and $n_{0}$ a reference number of
packets (for example the value at the lowest tested power). $i$ is the mode
index in the $j$-th group. In the case of a GRIN fiber, there are $2M=Q(Q+1)$
modes and polarizations, and $i=1,2,..M$; it is $i=1$ for $j=1$, $i=2,3$ for
$j=2$, $i=4,5,6$ for $j=3$ and so on. $g_{i}=[2,4,4,6,6,6,..,2Q]$,
$\epsilon_{i}=\beta_{i}-\beta_{j=Q}$ are the differential modal eigenvalues
with $\beta_{i}$ (m-1) the propagation constants. Eq. 1 was obtained under the
reasonable condition
$n_{0}>>\exp{\big{(}-\frac{\mu^{\prime}+\epsilon_{i}}{T}\big{)}}-1$; it is
able to reproduce the modal power distribution both in linear regime, where
RMC is mostly responsible for power exchange among modes, and in nonlinear
regime, where inter-modal four-wave mixing (IM-FWM) dominates. Chemical
potential $\mu^{\prime}$ (m-1), temperature $T$ (m-1) and the normalized power
$\gamma$ are three degrees of freedom for fitting Eq. 1 to the experimental
data, with the only constraint for $\gamma$ to scale with the input power, and
to respect the conservation law $\sum_{i=1}^{M}\lvert f_{i}\rvert^{2}=1$.

The thermodynamic validity of Eq. 1 can be stressed using the state equation
in Methods Sec. 6.1, Eq. 18, which can be reformulated as

$SE=\sum_{j=1}^{Q}\beta_{j}j\lvert
f_{j}\rvert^{2}+\mu^{\prime}-\beta_{j=Q}+(2M-Q)\frac{T}{\gamma}=0\mathchar
24635\relax\;$ (2)

an experimental error upon the state equation can be calculated as

$\epsilon_{SE}=\frac{SE}{\sum_{j=1}^{Q}\beta_{j}j\lvert
f_{j}\rvert^{2}-\Big{(}\mu^{\prime}-\beta_{j=Q}+(2M-Q)\frac{T}{\gamma}\Big{)}}.$
(3)

In Eq. 2, the first term in the sum (the normalized energy) must be constant;
$\mu^{\prime}$ and $T/\gamma$ must scale in order to maintain constant the
remaining part of the equation. Not necessarily the parameters $T/\gamma$ and
$\mu^{\prime}$ must remain separately constant.

In short fiber spans at low power, both RMC and IM-FWM are negligible; we can
assume equal modal distributions at the input and output of the fiber $\lvert
f_{j}^{(in)}\rvert^{2}=\lvert f_{j}^{(out)}\rvert^{2}$. Since the input modal
distribution is unchanged at high power, Eq. 2 can be used, together with
$\lvert f_{j}^{(out)}\rvert^{2}$ measured at low power, to predict the
thermodynamic parameters of the system at high power, when thermalization, or
condensation, is achieved pourbeyram2022direct . However, in long spans of
fiber at low power, RMC causes $\lvert f_{j}^{(in)}\rvert^{2}\neq\lvert
f_{j}^{(out)}\rvert^{2}$; thermodynamic parameters are found in this case by
fitting Eq. 1 to the experimental distribution at a given power; Eq. 2 is used
to verify the thermodynamic validity of the fit, but not for parameter
predictions from low power experiments.

The validity of Eq. 1 is also subject to the constraints of negligible
variations of the power $P$ and the energy $U$ (Sec. 6.1). Linear losses are
negligible, at the telecom wavelength, for fiber length up to few kilometers.
Raman nonlinearity is responsible for a frequency shift of the pulse spectrum,
and for the consequent variation in the packet energy; in soliton propagation,
Raman self-frequency shift (RSFS) Gordon1986  induces a red-shift of the
spectrum without distortion; assuming a tolerance of 5% on the energy change,
a red-shift of 80 nm at 1550 nm wavelength can be tolerated in the
experiments. Chromatic dispersion broadens pulses in the linear regime,
considerably reducing the pulse peak power $P$; however, in the quasi-soliton
and soliton regimes, the propagating train of pulses conserve their
pulsewidth, and negligible changes of power can be assumed over large portions
of the fiber. RMC is also responsible for negligible energy variations; let us
assume, by way of example, the propagation at 1550 nm in a GRIN fiber, and
that RMC is responsible for a complete power transfer from group 10 to the
fundamental mode (this is, of course, a worst case); the fractional energy
change equals the one of the propagation constants
$(\beta_{1}-\beta_{10})/\beta_{1}=0.0078$. In practical cases, strong RMC
induces negligible energy variations, comparable to the linear losses from few
meters of fiber. We can assume then, that the presence of strong RMC does not
invalidate the thermodynamic approach.

Experimental modal distributions that are suitably fitted by Eq. 1 correspond
to extrema of the entropy $S$ denoting the achievement of a steady state,
eventually characterized by a local or global condensation; power fluctuations
related to local condensation do not invalidate the equation, which is used to
fit the power distribution of the modal groups as a whole; in the next
section, Eq. 1 will be used for this purpose.

##  3 Experiments

In order to determine the mode power distribution at the output of GRIN MM
fibers, we used the mode-decomposition method introduced in Zitelli:23 .
Details about our experimental setup are provided in Methods Sec. (6.2). We
have progressively varied the input pulse energy $E_{in}$, so that we could
explore the whole spatio-temporal propagation regime from the linear to the
nonlinear case. We used 250 fs full-width-at-half-maximum (FWHM) pulses at
1400 nm, with 100 KHz repetition rate. The input beam was circularly
polarized, and coupled with a waist of 13 $\mu$m and a 10 $\mu$m lateral shift
with respect to the fiber axis, in order to enhance the population of higher-
order modes (HOM), and minimize power exchange between polarizations. We used
commercial OM4 GRIN (Thorlabs GIF50E) fibers, with lengths of 1 m, 830 m and 5
km, respectively, spooled on a support with radius of curvature larger than 8
cm.

Fig.2 shows the normalized instantaneous output power (left) and the near-
field (right) after 830 m of GRIN. Modal dispersion leads to the observed time
delay among pulses carried by different groups of degenerate modes. Note that
the sub-pulses are equally spaced in time, owing to the equal spacing of the
mode propagation constants in GRIN fibers. The pulse carried by the HOMs has
the largest delay, hence it appears in the trailing tail of the output
waveform.

Using a fiber which is long enough to temporally separate from different mode
groups permits to directly measure the output mode power distribution. This
results from the combined action of linear and nonlinear mode coupling, which
occurs before pulse splitting. The linear and nonlinear interaction among
different modal groups is mostly concentrated in the first portion of the
fiber, where pulses overlap in time. Hence, negligible linear losses can be
assumed over the interaction length, thus permitting a correct thermodynamic
interpretation of the results.

In the linear regime, we estimate that chromatic dispersion (CD) broadens the
individual pulses in each group up to a 110 ps duration after 830 m (when
considering the nominal -12 ps2/km chromatic dispersion). Modal dispersion
produces a time delay of 206 ps between adjacent group pulses. Optical pulse
energy is varied from linear regime (0.1 nJ) to soliton regime (4.5 nJ).
Quasi-soliton regime is obtained from 1.0 to 4.1 nJ energy, characterized by
trains of short pulses. Raman soliton regime is reached at 4.5 nJ: no Raman
delay or spectral shift was observed up to 4.1 nJ; hence, a conservative
system is assumed up to 4.1 nJ.

In the linear regime (Fig.2a), the output distribution of energy among the
mode groups is determined by linear disorder. For $E_{in}=0.81$ nJ, the
effective length of the fiber (determined by taking into account both the weak
linear loss and the rapid pulse peak power decrease owing to chromatic
dispersion-induced temporal broadening) is as short as $L_{eff}=12$ m. This
means that, for most of the distance, power exchange among pulses carried in
different modes is only determined by RMC.

In Fig.2a, c, e, we sample the photodiode traces at points corresponding to
the relative group delay of each modal group (see the orange circles), which
permits a mode power decomposition as described in Zitelli:23 . In this
approach, and considering the first $Q=10$ modal groups, corresponding to
$M=Q(Q+1)/2=55$ modes per polarization, one directly measures the mode group
powers $P_{j}$. The corresponding mean modal power fraction carried by each
mode in the group is then $\lvert f_{i}\rvert^{2}=2P_{j}/(g_{i}P_{tot})$, by
supposing power equipartition within each group (see Methods Sec 6.3).

Figure 2: Normalized group power fractions (left) with their corresponding
near-field intensities (right) after 830 m of GRIN fiber, for input pulse
energy of a-b) 0.81 nJ, c-d) 2.45 nJ, e-f) 4.10 nJ.

Fig.2a shows that, in the linear propagation regime (i.e., for $E_{in}\leq
0.81$ nJ) most of the pulse energy is carried by the first 6 mode groups under
the experimental coupling conditions, and mostly in the second group. On the
other hand, we can see in Fig.2c that for $E_{in}=2.45$ nJ, the lower-order
modes attract power from the HOMs as a consequence of the increasing IM-FWM.
At this energy levels, a strongly nonlinear quasi-soliton regime occurs.

At 4.10 nJ pulse energy (Fig.2e), the fundamental mode has attracted half of
the total power, and the propagation is approaching to a multimode soliton.
Still, the soliton has not fully formed and does not show the characteristic
Raman delay of a full soliton regime Zitelli2021a ; for this reason, we may
neglect the presence of Raman scattering and other dissipative effects up to
about this level of input pulse energy.

In the experiments, we also observed that for $E_{in}\geq 4.6$ nJ (not shown),
a soliton is formed in the fundamental mode: for larger values of $E_{in}$,
the soliton is further delayed in time and splits away from the remaining
pulse carried by HOMs, because of the Raman soliton self-frequency shift
Gordon1986 ; Renninger2013 ; Wright:15 ; Zitelli_9887813 .

As shown in Figs.2b, d, the output mode power distribution leads to a
relatively broad and speckled output beam intensity profile at both low and
intermediate values of $E_{in}$. Whereas approaching to the soliton regime,
most of the energy is attracted by the fundamental mode; as a result, Fig.2f
shows that the beam brightness is substantially enhanced at its center: a
bell-shaped beam is formed with a waist close to that of the fundamental mode,
sitting on a background of HOM, a typical signature of beam self-cleaning
krupa2017spatial .

Figure 3: Average mode power fraction $\lvert f_{i}\rvert^{2}$ vs. modal
eigenvalue $\epsilon_{i}$, for input pulse energies $E_{in}$ equal to a) 0.16
nJ,(b) 0.82 nJ, (c) 2.45 nJ, and d) 4.10 nJ, respectively. Corresponding
analytic fits with the weighted BE distribution are shown by a solid orange
curve.

Thanks to the accuracy offered by the used decomposition method, Fig.3 reports
in log scale the average output power fractions $\lvert f_{i}\rvert^{2}$ from
Fig.2, versus their respective mode eigenvalues $\epsilon_{i}$. In Fig.3, the
input pulse energy ranges from $E_{in}=0.2$ nJ, where mode mixing originates
from RMC, up to $E_{in}=4.1$ nJ, where IM-FWM is the dominant mechanism for
transferring energy among nondegenerate modes krupa2017spatial .

Output modal distributions are fitted using Eq. 1, taking care to respect the
constraint for the normalized power $\gamma$. At low pulse energy (Fig.3a),
the weighted BE approximates to a straight line, in good agreement with
predictions obtained by numerically solving the so-called power-flow equations
Savovi2019PowerFI  Gloge:6774107  (see Sec 4.1). For pulse energy of 0.81 nJ
(Fig.3b), propagation is initially nonlinear, and later dominated by RMC. The
weighted BE properly fits the experimental data up to the 9-th group order,
for $\mu^{\prime}=-58.2$ mm-1, $T=8.71$ mm-1 and $\gamma=3.60$. In the figure,
we observe that group 2 has larger energy fraction respect to the fitting
equation. It will be shown later that a local condensation of energy into
lower groups is obtained at this power level.

On the other hand, Figs.3c, d show that, as soon as $E_{in}\geq 2.45$ nJ, the
population of the fundamental mode grows larger, so that it preferentially
acquires power from HOMs as a consequence of IM-FWM. In both cases, the
weighted BE fits the experimental data up to the soliton regime, obtaining
$\mu^{\prime}=-53.1(-52.2)$ mm-1, $T=8.71(10.14)$ mm-1 and $\gamma=10.8(18.0)$
for $E_{in}=2.45(4.10)$ nJ, respectively; the BE distribution denotes a new
equilibrium state that is induced by the strong nonlinearity, which is
cumulated by the pulse over the entire fiber length. Note that, for
$E_{in}=4.1$ nJ one has that $(\mu^{\prime}+\epsilon_{1})/T\simeq-0.13$, which
means that the RJ approximation to the BE law (i.e., $\lvert
f_{i}\rvert^{2}\propto-T/(\mu^{\prime}+\epsilon_{i})$) is appropriate around
the fundamental mode with $i=1$.

Figure 4: a) Experimental group power fraction $j\cdot\lvert f_{j}\rvert^{2}$
vs. group index $j$ for energies ranging between 0.16 nJ and 4.61 nJ. b)
Simulated group power fraction, in similar conditions to the experiments.

Table 1 shows the values of the optical temperature $T$, differential
potential $\mu^{\prime}$, factor $\gamma$ and state equation error
$\epsilon_{SE}$ for the experiment of Fig.3. For growing input energy
$E_{in}$, temperature $T$ increases from 7.11 mm-1 to 10.36 mm-1;
correspondingly, $\mu^{\prime}$ reduces from -65.42 mm-1 to -52.08 mm-1;
$\gamma$ scales proportionally to the input energy, and the error upon the
state equation decreases from 7.9% and 3.2% in the linear regime down to 0.3%
in the nonlinear regime.

Table 1: Thermodynamic parameters obtained fitting the data of Fig.3 using the weighted BE law. $E_{in}$ [nJ] |  $T$ [mm-1] |  $\mu^{\prime}$ [mm-1] | $\gamma$ | $\epsilon_{SE}$
---|---|---|---|---
0.16 | 7.11 | -65.42 | 0.72 | 0.0793
0.41 | 7.76 | -60.92 | 1.80 | 0.0318
0.82 | 8.71 | -58.19 | 3.60 | 0.0177
1.23 | 9.32 | -56.52 | 5.40 | 0.0122
1.64 | 9.05 | -54.74 | 7.20 | 0.0085
2.06 | 9.53 | -54.17 | 9.00 | 0.0069
2.45 | 8.71 | -53.10 | 10.80 | 0.0058
3.29 | 9.75 | -52.81 | 14.40 | 0.0041
4.11 | 10.14 | -52.29 | 18.00 | 0.0035
4.61 | 10.36 | -52.08 | 20.16 | 0.0035
|  |  |  |

A detailed information on how the relative power fraction for groups of
nondegenerate modes, or $j\cdot\lvert f_{j}\rvert^{2}$, evolves as a function
of group index $j$ for several $E_{in}$ is given in Fig.4a. In the figure, it
appears that in the intermediate nonlinear regime (0.82 nJ), group 2 locally
attract energy at the expense of the HOMs; at 1.65 nJ and 3.29 nJ, local
condensation into group 4 is observed. When approaching to the soliton regime
(4.61 nJ), the fundamental mode power grows up larger at the expenses of
groups 2-4 and the HOMs, until it carries more than 50% of the total output
pulse energy. For $E_{in}>4.6$ nJ, a Raman soliton or a train of solitons
form, which can be spectrally isolated from the residual pulses; in this
regime (not shown), experiments confirm that more than 90% of the pulse energy
is concentrated to the ground state.

Fig.4b reports similar results from the numerical simulations explained in
section 4.2. Here, it is evident the local condensation into group 2 at 1.0
nJ, and into group 3 at 3.0 nJ, followed by a global condensation into group 1
at 5.0 nJ.

In the supplementary Sec. A, a validation of the weigthed BE law is reported
against independent experimental data from pourbeyram2022direct .

##  4 Simulations

###  4.1 Linear disorder

In the linear propagation regime, modal power exchanges originate from RMC. A
well-know model to simulate linear coupling among modes in a GRIN fiber is
provided by the power-flow equations Savovi2019PowerFI  Gloge:6774107  (see
Sec. 6.3). The model describes a bi-directional flow of power, from group $j$
down to group $j-1$ and up to group $j+1$. After a number of consecutive
integration steps, a cascading effect is produced, leading to power flow into
non-adjacent groups. As discussed in Methods Sec. 6.3, the mode coupling
coefficients are not symmetrical in the two directions, resulting into a
preferential transfer of power from HOMs into low-order modes of the fiber.

An example of numerical solution of the power-flow equations is provided in
Fig.5. Here we simulated the propagation of the first 10 mode groups at 1400
nm, over 830 m of GRIN fiber. The linear loss is $\alpha_{0}=2.6$ dB/km, and
the coupling coefficient is $D=0.043$ m-1; we also included the presence of
weak modal losses via the coefficient $A=2\times 10^{-5}$ m-1. The modal
content at the fiber input is assumed to be equi-distributed; hence, the group
power increases linearly with group index $j$.

Fig.5a shows the evolution of the power of the groups vs. distance: a net flow
of power in the direction of the lower-order groups is observed. As a result,
the output mean modal power fraction $\lvert f_{i}\rvert^{2}$ vs. the mode
eigenvalue progressively increases (Fig.5b). This distribution can be fitted
well by the weighted BE law with $\mu^{\prime}=-54.6$ mm-1, $T=31.1$ mm-1 and
$\gamma=70.5$, which describes the RMC-induced steady-state. Figure also shows
that a simple exponential function, of the type $\lvert
f_{i}\rvert^{2}=a\exp{(-b\epsilon_{i})}$ cannot properly fit the HOM groups
for low losses.

Figure 5: a) Power of the mode groups vs. distance, resulting from a power-
flow simulation in the linear regime. b) Corresponding input and output mean
modal power fractions vs. the modal eigenvalues.

###  4.2 Linear disorder and nonlinearity

In order to provide additional insight into the mode coupling dynamics leading
to our experimental results, it is necessary to introduce a model which takes
into account the presence of both linear and nonlinear mode coupling. To this
end, we carried out extensive numerical simulations by using the coupled-mode
GNLSEs Poletti2008a , including wavelength-dependent linear and modal losses,
and a new model for RMC which is derived from the power-flow equations of a
GRIN fiber Savovi2019PowerFI , as explained in the Method Sec. (6.4). To save
computation time, only $M=28$ modes were propagated over 100 m of GRIN fiber;
we considered an input 250 fs pulse at 1400 nm, with the same coupling
conditions as in the experiments. Modes 1 to 28 correspond to the Laguerre-
Gauss modes $LG_{01}$, $LG_{11e}$, $LG_{11o}$, $LG_{21e}$, $LG_{21o}$,…, and
$LG_{04}$, respectively. We set the RMC coupling coefficient $D=0.003$ m-1,
which introduces a significant amount of linear disorder over the considered
distance. The RMC step was equal to $L_{c}=6$ mm, to ensure appropriate
simulation accuracy. The fiber parameters (second and third-order chromatic
dispersion, modal dispersion, Kerr and Raman nonlinearity, linear losses)
correspond to typical values for an OM4 GRIN fiber at 1400 nm.

Figure 6: Simulated output modal power for $E_{in}$ equal to a) 0.02 nJ, b)
3.0 nJ and c) 5.0 nJ, respectively, with the same input coupling conditions as
in the experiment.

Fig 6a shows the simulated temporal dependence of the output modal power in
the linear regime ($E_{in}=0.02$ nJ). As can be seen, in each mode group the
input pulses are broadened up to 13 ps by chromatic dispersion. Moreover,
pulses in different mode groups separate in time by the action of inter-modal
dispersion. On the other hand, Fig.6b shows that in the nonlinear propagation
regime ($E_{in}=3$ nJ), pulses within each mode group are significantly
compressed in time by the combined action of SPM and anomalous dispersion;
pulses do not experience Raman delay, meaning that a full soliton regime is
not achieved yet. IM-FWM transfers power principally to the lower-order modes.

Finally, Fig.6c shows that, in the soliton regime ($E_{in}=5$ nJ), ultrashort
(170 fs duration) Raman delayed soliton pulses are formed; nearly 90% of the
power is attracted in the fundamental mode of all travelling pulses. A train
of fundamental solitons results after long distance.

The resulting numerically simulated mode group power fraction $j\cdot\lvert
f_{j}\rvert^{2}$ distribution at the GRIN fiber output is shown in Fig.4b for
the first 4 mode groups. As can be seen, in the linear propagation regime
($E_{in}=0.02$ nJ or 1.0 nJ), group power drops following a convex curve (in
log scale). Whereas whenever Raman solitons are formed ($E_{in}=5.0$ nJ), the
group power decreases following a concave curve: this is associated with a
nonlinear irreversible energy condensation into the fundamental mode
Zitelli2021 . At intermediate energy (2.0 to 3.0 nJ), one would expect a
linear decreasing shape of the group power curve; to the contrary, local
energy condensations are observed into groups 2 and 3.

It must be said that the weighted BE is able to fit the output modal
distributions (not shown) obtained from the simulations of Fig.4b, denoting
the achievement of steady states.

##  5 Discussion and Conclusions

The mean modal power fractions $\lvert f_{i}\rvert^{2}$ of Fig.3 show the
achievement of steady states in the quasi-soliton regimes, because the
weighted BE properly fits the distributions reached after long distances, as
discussed in Sec. 6.1. Power-flow simulations in the linear regime, after long
propagation distances, provide distributions properly fitted by the weighted
BE as well, confirming that Eq. 1 is a reliable method to analyze the
achievement of steady states even at low power.

By representing the mean modal power fraction from the same modal group,
against the input pulse energy (Fig.4a), the fundamental mode (group 1)
reaches the highest fraction in the soliton regime, where IM-FWM combined with
RMC is responsible for power redistribution. For intermediate energy levels,
between 0.4 and 2.0 nJ, groups 2 and 3 increase their fraction at the expense
of the HOMs. A local attraction of energy to the lower-order groups is
observed.

The experimental findings of Fig.4a are explained by observing results of the
numerical simulations in Figs.6 and 4b. In all regimes (low power, quasi-
soliton and soliton), RMC transfers power to modes within the same or of the
neighbouring groups; in the linear regime, this process has the macroscopic
effect to produce a BE distribution. However, in the quasi-soliton and soliton
regimes, the Kerr nonlinearity and anomalous chromatic dispersion considerably
shorten pulses; RMC diffuses power to the lower-order modes, a process which
is boosted by the IM-FWM.

Globally, thanks to the interplay of RMC and IM-FWM, energy clusters are
formed in the lower-order groups in the quasi-soliton regime (Figs.6b, 4 and
1b); a train of quasi-soliton pulses is produced, each composed by the modes
of one group $j$, plus a fraction of modes belonging to the lower adjacent
groups and of the fundamental mode; this results in the promotion of the first
3-4 modal groups at the output, and power clustering. In the framework of such
interplay, RMC invalidates the achievement of a global condensation state as
described by a RJ law PhysRevLett.129.063901 ; however, the formation of
steady states associated to a local condensation is still possible, a process
which is characterized by a weighted BE modal distribution.

When increasing the pulse energy up to the soliton regime (Figs.6c and 1c),
global modal condensation to the ground state is observed in all splitted
pulses; nearly 90% of the power is attracted to the fundamental mode, as a
consequence of the interplay of RMC and IM-FWM. Pulsewidth reduces to 170 fs,
which is typical of propagated walk-off solitons Zitelli2021a . A train of
fundamental solitons is eventually produced by a peculiar fission mechanism
mediated by modal dispersion; then solitons are affected by Raman soliton
self-frequency shift Zitelli_9887813 .

In thermodynamic terms, experimental and numerical data confirm that the
linear regime distribution describes a gas of energy packets, which evolves
into a local condensed, or glassy state in the intermediate nonlinear regime,
and then into a globally condensed solid state in the soliton regime
Jung2019MeasuringTO . All of the observed distributions correspond to steady
states, being their modal distributions properly fitted by the weighted BE.

##  6 Methods

###  6.1 Theory

In the main text it was considered an optical multimode system including $Q$
groups of degenerate modes, distributed over $g_{j}/2$ modes and 2
polarizations, with $j=1,2,..,Q$; hence, $g_{j}$ is the degeneracy over two
polarizations; in the special case of a GRIN fiber, it is $g_{j}=2,4,6,..,2Q$;
the number of modes and polarizations is $2M=Q(Q+1)$. Following the procedure
proposed in Wu2019 , the population $n_{j}$ of energy packets into the j-th
group over 2 polarizations, leads to a total number of combinations among
groups and polarizations

$W=\prod_{j=1}^{Q}\frac{(n_{j}+g_{j}-1)!}{(g_{j}-1)!n_{j}!}.$ (4)

System entropy is defined as $S=ln(W)$; by applying the Stirling
approximation, valid for $n_{j}+g_{j}-1>>1$ we obtain

$S=\sum_{j=1}^{Q}(n_{j}+g_{j}-1)\big{[}ln(n_{j}+g_{j}-1)-1\big{]}-(g_{j}-1)\big{[}ln(g_{j}-1)-1\big{]}-n_{j}\big{[}ln(n_{j})-1\big{]}\mathchar
24635\relax\;$ (5)

Global system thermalization corresponds to the extremization of Eq. 5 using
Lagrange multipliers for the total number of particles $N=\sum_{j}n_{j}$ and
total normalized energy $E=\sum_{j}\beta_{j}n_{j}$, being $\beta_{j}$ the
modal propagation constants supposed equal for the degenerate modes

$\frac{\partial}{\partial
n_{l}}\Big{[}S+\sum_{j=1}^{Q}\Big{(}a^{\prime}n_{j}+b^{\prime}\beta_{j}n_{j}\Big{)}\Big{]}=0,$
(6)

which provides, with no approximations

$ln\Big{(}1+\frac{g_{j}-1}{n_{j}}\Big{)}+a^{\prime}+b^{\prime}\beta_{j}=0\mathchar
24635\relax\;$ (7)

Eq. 7 is valid for

$\frac{n_{j}}{g_{j}-1}=\frac{1}{\exp{[-(a^{\prime}+b^{\prime}\beta_{j})]}-1}.$
(8)

We can choose $a^{\prime}=\mu/(Tn_{0})$ and $b^{\prime}=1/(Tn_{0})$, with
$n_{0}$ a reference number of energy packets (for example the value at the
lowest tested power) and $N=\gamma n_{0}$. $T$ (1/m) is an optical temperature
and $\mu$ (1/m) a chemical potential. We these choices, Eq. 8 can be written
as

$n_{j}=\frac{(g_{j}-1)}{\exp{\Big{(}-\frac{\mu+\beta_{j}}{Tn_{0}}\Big{)}}-1}.$
(9)

An alternative development of Eq. 8 consists in replacing $a^{\prime}$ and
$b^{\prime}$ with non-factorizable constants $a$ and $b$ defined as

$-(a+b\beta_{j})=\ln{\Big{[}\frac{1}{n_{0}}\exp{\Big{(}-\frac{\mu+\beta_{j}}{T}\Big{)}}-\frac{1}{n_{0}}+1\Big{]}}\simeq\frac{1}{n_{0}}\Big{[}\exp{\Big{(}-\frac{\mu+\beta_{j}}{T}\Big{)}}-1\Big{]}.$
(10)

The approximation in Eq. 10 is valid for $\lvert a+b\beta_{j}\rvert<<1$, which
is less than 10-5 in the experiments. By replacing into Eq. 8 we obtain

$n_{j}=\frac{n_{0}(g_{j}-1)}{\exp{\Big{(}-\frac{\mu+\beta_{j}}{T}\Big{)}}-1}.$
(11)

The particular choice of $a$, $b$, $a^{\prime}$ and $b^{\prime}$ conserves the
system’s energy and power; in fact it results, for $a^{\prime}$, $b^{\prime}$

$\sum_{j=1}^{Q}a^{\prime}n_{j}=-\sum_{j=1}^{Q}\frac{\mu(g_{j}-1)}{\mu+\beta_{j}},$
(12)

$\sum_{j=1}^{Q}b^{\prime}\beta_{j}n_{j}=-\sum_{j=1}^{Q}\frac{\beta_{j}(g_{j}-1)}{\mu+\beta_{j}}\mathchar
24635\relax\;$ (13)

and

$\sum_{j=1}^{Q}a^{\prime}n_{j}+b^{\prime}\beta_{j}n_{j}=-\sum_{j=1}^{Q}(g_{j}-1)=Q-2M.$
(14)

For $a$, $b$ it also results

$\sum_{j=1}^{Q}(a+b\beta_{j})n_{j}=-\sum_{j=1}^{Q}(g_{j}-1)=Q-2M,$ (15)

hence, the choices of $a$, $b$ or $a^{\prime}$, $b^{\prime}$ are equivalent in
terms of power and energy conservation.

When passing from Eq. 9 to Eq. 11 by replacing the $a^{\prime}$, $b^{\prime}$
with $a$ and $b$, the extremization of the entropy is performed using non-
factorizable multipliers directly related to the modal group eigenvalues. The
entropy extrema assume a local significance, because both $\mu$ and $T$ are
now related to the particular set of $\beta_{j}$ and cannot be isolated. The
total system’s power and energy is conserved in the same way as for the global
procedure of Eq. 9. The solution Eq. 11 is able to indicate the achievement of
energy minima by the individual modal groups, such as locally condensed
states; the equation is suitable to fit experimental modal distributions as a
whole, regardless of the power fluctuations caused by local condensates.

We now pass to the $i$-th mode in the $j$-th group; in the case of a GRIN
fiber, there are $2M=Q(Q+1)$ modes and polarizations, and $i=1,2,..M$; it is
$i=1$ for $j=1$, $i=2,3$ for $j=2$, $i=4,5,6$ for $j=3$ and so on, and
$g_{i}=[2,4,4,6,6,6,..,2Q]$. The mean modal power fraction over two
polarizations is obtained as $\lvert f_{i}\rvert^{2}=2n_{j}/(\gamma
n_{0}g_{i})$, with $\gamma=N/n_{0}$; $\gamma>1$ must scale with the
experimental power.

Introducing the differential eigenvalues $\epsilon_{i}=\beta_{i}-\beta_{j=Q}$,
referred to the higher mode, and $\mu^{\prime}=\mu+\beta_{j=Q}$, eq. 11
provides the weighted Bose-Einstein (BE) modal distribution

$\lvert
f_{i}\rvert^{2}=\frac{2(g_{i}-1)}{g_{i}\gamma}\frac{1}{\exp\big{(}-\frac{\mu^{\prime}+\epsilon_{i}}{T}\big{)}-1}.$
(16)

The constraints for $\gamma$ are to scale with the input power, and to respect
the conservation law $\sum_{i=1}^{M}\lvert f_{i}\rvert^{2}=1$.

Let $P$ and $P_{0}$ be the optical powers corresponding to $N$ and $n_{0}$
energy packets (for an optical pulse, it is the peak power or the energy);
hence, $\gamma=N/n_{0}=P/P_{0}$. The system’s internal energy is
$U=-\sum_{j}\beta_{j}n_{j}P/N$ (W/m); the power $P=\sum_{j}n_{j}P/N$ (W). From
Eqs. 14 and 15 we get

$\sum_{j=1}^{Q}(a+b\beta_{j})n_{j}=Q-2M=\frac{\mu
N}{Tn_{0}}+\frac{1}{Tn_{0}}\Big{(}-\frac{UN}{P}\Big{)},$ (17)

which provides the state equation

$U-\mu P=(2M-Q)P_{0}T.$ (18)

From the local extremization problem, providing Eq. 11, it is not possible to
derive the well-known Rayleigh-Jeans (RJ) distribution Wu2019
podivilov2022thermalization  under reasonable approximations. However, from
the global extremization solution, Eq. 8 with the definitions of $a^{\prime}$
and $b^{\prime}$, we obtain

$\lvert f_{i}\rvert^{2}=\frac{2(g_{i}-1)}{g_{i}\gamma
n_{0}}\frac{1}{\exp\big{(}-\frac{\mu^{\prime}+\epsilon_{i}}{Tn_{0}}\big{)}-1}\mathchar
24635\relax\;$ (19)

the presence of $N=\gamma n_{0}$ at the denominator of Eq. 19 makes the
equation suitable for fitting the experimental data only for $N<10$; for this
reason, it is unusable. The problem is removed under the assumption
$\lvert\mu^{\prime}+\epsilon_{i}\rvert<<\lvert Tn0\rvert$, which brings
directly to the RJ distribution

$\lvert
f_{i}\rvert^{2}=-\frac{2(g_{i}-1)}{g_{i}\gamma}\frac{T}{\mu^{\prime}+\epsilon_{i}}.$
(20)

In terms of fractional power $\lvert c_{i}\rvert^{2}=P\lvert f_{i}\rvert^{2}$
in [W], by choosing $T^{\prime}=P_{0}T$ (W/m), Eq. 20 provides.

$\lvert
c_{i}\rvert^{2}=-\frac{2(g_{i}-1)}{g_{i}}\frac{T^{\prime}}{\mu^{\prime}+\epsilon_{i}}.$
(21)

The RJ law demonstrates to be a good choice for describing experiments
characterized by global condensation states, such as the well-known self-
cleaning experiments krupa2017spatial ; mangini2022statistical ;
podivilov2022thermalization , while the more general weighted BE, Eq. 16, is
suitable to fit the experimental data also in cases where local condensed
states are achieved.

###  6.2 Experimental setup

Optical pulses at 1400 nm wavelength and 250 fs pulsewidth are generated by an
optical parametric amplifier (OPA) fed by a femtosecond Yb laser, at 100 kHz
repetition rate, and using a narrow-band optical filter. The input beam is
attenuated, linearly polarized and passed through a $\lambda/4$ waveplate, in
order to generate a circular state of polarization. The Gaussian shaped beam
is then injected into variable length spans of OM4 GRIN fiber (1 m, 830 m and
5 km), with a waist $w_{0}$ of approximately 13 $\mu$m. The corresponding
self-imaging induced beam compression factor
$C=2z_{p}/(\pi\beta_{0}w_{0}^{2})=0.305$, where $\beta_{0}=2\pi
n_{0}/\lambda$, $z_{p}=0.55$ mm the self-imaging period, and $n_{0}=1.46$ is
the core index Karlsson:92 ; Ahsan:19 ; the effective beam waist for
determining the nonlinear coefficient is $w_{e}=\sqrt{C}w_{0}=7.27\mu$m
Ahsan:18 ; Ahsan:19 , close to the fundamental mode waist. The circular state
of polarization is used at the input in order to minimize power exchanges
among polarizations. The input beam is laterally shifted with respect to the
fiber axis by 10 $\mu$m, in order to increase the contribution from higher-
order modes.

Linear losses were measured as $\alpha=6.0\text{x}10^{-4}$ m-1. According to
Gloge theory Gloge:72 , bending losses remain negligible up to the first 10
mode groups ($\alpha_{j}<=8\text{x}10^{-10}$ dB/km for group $j=1,2,..,10$,
and $\alpha>>30$ dB/km for group $j>=11$); in agreement with the theory, 10
modal groups could be observed at the output of both 830 m and 5 km fiber
spans.

Modal dispersion is responsible for the time delay of the different groups;
measured delay among groups is 206 ps over 830 m, and 850 ps over 5 km. Hence,
mode groups are time-resolved at the output and easily measurable after 830 m,
once the modes have interacted linearly and nonlinearly along the fiber. Delay
among groups was found to vary with distance as $\Delta t_{j}=1.02z^{0.79}$;
RMC affecting the experiment is intermediate between weak (where modal delay
scales with $z$) and strong regime (where it changes with $\sqrt{z}$)
Marcuse1973LossesAI .

At the fiber output, near-field is imaged on an InGaAs camera (Hamamatsu
C12741-03); the beam is also directed to a real-time multiple octave spectrum
analyzer with a spectral detection range of 1100-5000 nm (Fastlite Mozza). The
output pulse instantaneous power is detected by a fast photodiode (Alphalas
UPD-35-IR2-D) and a real-time oscilloscope (Teledyne Lecroy WavePro 804HD)
with 30 ps overall response time. An intensity autocorrelator (APE pulseCheck
50) with femtosecond resolution is also used for the characterization of the
input pulses. Input and output power are measured by a power meter with $\mu$W
resolution.

Traditional 2D modal decomposition methods Flamm2012  are not suitable for the
analysis of the output near-field after hundreds of meters of pulse
propagation, because they do not account for: (i) the phase chirp which is
induced by chromatic dispersion of pulses carried by the different modes; (ii)
the phase delay among modes due to the modal dispersion; (iii) laser-induced
phase noise; (iv) random phase differences among modes, introduced by the RMC.
In this work, we use the 3D modal reconstruction method proposed in Zitelli:23
; the mode group power is measured from the instantaneous power detected by a
fast photodiode; comparison to the measured near-field is performed after 3D
reconstruction. Such method has provided an accurate estimate of the modal
distributions reached in long spans of GRIN, both in linear and nonlinear
regime, up to the 10-th modal group (55 modes per polarization).

###  6.3 Power-flow numerical model

In the linear regime, RMC is properly modelled by the well-known power-flow
diffusive equations. If $P_{j}$ is the power of the $j$-th mode group, the
power exchange among adjacent modal groups is described by Savovi2019PowerFI

$P_{j}(z+L_{c})=\big{(}\frac{DL_{c}}{\Delta m^{2}}-\frac{DL_{c}}{2m\Delta
m}\Big{)}P_{j-1}(z)+\\\ \Big{(}1-\frac{2DL_{c}}{\Delta
m^{2}}-(\alpha_{0}+Am^{2})L_{c}\Big{)}P_{j}(z)+\big{(}\frac{DL_{c}}{\Delta
m^{2}}+\frac{DL_{c}}{2m\Delta m}\Big{)}P_{j+1}(z),$ (22)

with $m(j)=j-1$, $D$ and $\alpha_{0}$ (1/m) the coupling coefficient and
linear losses, respectively, $A$ the modal loss coefficient, $\Delta m=1$ the
modal step and $L_{c}$ the RMC integration step. In Eq. 22, power flows in
both directions, from group $j$ down to group $j-1$ and up to $j+1$; after
consecutive integration steps, a cascading effect is produced, causing the
coupling of non-adjacent groups.

Mode coupling into groups is neglected, because it is so fast that statistical
modal equipartition into groups can be assumed; it is then meaningful to
calculate the mean modal content into groups as $\lvert
f_{i}\rvert^{2}=2P_{j}/(g_{i}P_{tot})$.

####  6.3.1 Modal power equipartition

It is possible to demonstrate, using the alternative Gibb’s definition of
entropy doi:10.1119/1.1971557 , that modal equipartition does not generally
apply among different groups. Given, at steady state, $p_{i}=<\lvert
f_{i}\rvert^{2}>$ the modal occupation probability, and
$\lambda_{1},\lambda_{2}$ two Lagrange multipliers, extremization of the
entropy, while power $P$ and energy $E$ are conserved, reads

$\frac{\partial}{\partial p_{k}}\Big{[}-\sum_{i=1}^{M}p_{i}\log
p_{i}+\lambda_{1}\Big{(}\sum_{i=1}^{M}p_{i}-1\Big{)}+\lambda_{2}\Big{(}\sum_{i=1}^{M}p_{i}\beta_{i}-E/P\Big{)}\Big{]}=0\mathchar
24635\relax\;$ (23)

from Eq. 23, if only the power is conserved ($\lambda_{2}=0$), we obtain

$\sum_{i=1}^{M}p_{i}=\sum_{i=1}^{M}\exp{(\lambda_{1}-1)}=M\exp{(\lambda_{1}-1)}=1,$
(24)

which brings to the modal equipartition $p_{i}=1/M$. However, when considering
also the conservation of the energy, Eq. 23 provides

$\sum_{i=1}^{M}p_{i}=\sum_{i=1}^{M}\exp{(\lambda_{1}+\lambda_{2}\beta_{i}-1)}=1,$
(25)

which provides equipartition only for degenerate modes, with same $\beta_{i}$.

###  6.4 GNLSE numerical model

In the nonlinear (and linear) regime, numerical simulations use the coupled
GNLSEs Poletti:08 , modified to include modal and wavelength-dependent losses,
and linear random-mode coupling (RMC) for mode $p$

$\frac{\partial A_{p}(z,t)}{\partial
z}=i(\beta_{0}^{(p)}-\beta_{0})A_{p}-(\beta_{1}^{(p)}-\beta_{1})\frac{\partial
A_{p}}{\partial
t}+i\sum_{n=2}^{4}\frac{\beta_{n}^{(p)}}{n!}\big{(}i\frac{\partial}{\partial
t}\big{)}^{n}A_{p}-\frac{\alpha_{p}(\lambda)}{2}A_{p}+\\\
+i\sum_{m}q_{mp}A_{m}+in_{2}k_{0}\sum_{l,m,n}S_{plmn}\big{[}(1-f_{R})A_{l}A_{m}A_{n}^{*}+f_{R}A_{l}[h_{R}*(A_{m}A_{n}^{*})]\big{]}.$
(26)

In Eq. 26, $\beta_{n}^{(p)}$ is the n-th order dispersion term (modal and
chromatic) for mode $p$, $\alpha_{p}(\lambda)$ the modal and wavelength-
dependent loss coefficient, $n_{2}$ (m2/W) the nonlinear index coefficient
multiplying the Kerr and Raman terms, $S_{plmn}$ is an overlap integral among
modes, accounting for IM-FWM, and $q_{mp}$ is the linear RMC coupling
coefficient, from mode $m$ to $p$, coming from the power-flow equations model
Savovi2019PowerFI

$q_{mp}=\left\\{\begin{array}[]{cc}\big{[}\frac{DL_{c}}{p-1}\big{(}1-\frac{1}{2(p-1)}\big{)}\big{]}^{1/2}&\text{from
modes }m\text{ with }g_{m}=g_{p}-1\\\
\big{[}\frac{DL_{c}}{p+1}\big{(}1+\frac{1}{2(p-1)}\big{)}\big{]}^{1/2}&\text{from
modes }m\text{ with }g_{m}=g_{p}+1,\end{array}\right.$ (27)

being $g_{p}$ the degeneracy of modes, $p=1,2,..,M$, $D$ (m-1) the RMC
coupling coefficient, and $L_{c}$ the RMC numerical step. Degenerate modes are
not accounted for in Eq. 27, because their coupling is so fast that power
equipartition can be assumed into groups.

In the main text, fiber parameters at $\lambda=1400$ nm are: $\beta_{2}=-11.8$
ps2/km, $\beta_{3}=0.102$ ps3/km for the fundamental mode, $\alpha=2.6$ dB/km,
negligible modal loss, $n_{2}=2.7\text{x}10^{-20}$ m2/W, $D=0.003$ (1/m),
$L_{c}=6$ mm, $f_{R}=0.18$, $S_{plmn}$ calculated from modal overlap
integrals.

———————

Acknowledgments

The authors wish to thank C. Conti and G. D’Aguanno for valuable discussions
about glass states, S. Savovic’ for information regarding the power-flow
model, M. Ferraro, D. Christodoulides, G. Pyrialakos and G. Steinmeyer for
discussion regarding optical thermodynamics, and F. Wise and L. Wright for
making freely available the open-source parallel numerical mode solver for the
coupled-mode nonlinear Schrödinger equations Wright2018a , and the test data
used in supplementary Sec. A; a modified version of the software was used for
the numerical simulations.

## Declarations

  * •

Funding

European Union under the Italian National Recovery and Resilience Plan (NRRP),
Rome Technopole Flagship Project 5.

HORIZON EUROPE European Research Council (101081871).

NextGenerationEU, partnership on “Telecommunications of the Future”
(PE00000001 - program “RESTART”).

  * •

Conflict of interest/Competing interests

The authors declare no conflict of interest.

  * •

Ethics approval Not applicable

  * •

Consent to participate Not applicable

  * •

Consent for publication Not applicable

  * •

Availability of data and materials

Data available under request.

  * •

Code availability Not applicable

##  Appendix A Validation of the weighted BE law against independent
experiments

In order to provide an independent validation of the weighted BE law
introduced in Sec. 2, we referred to the tests performed in
pourbeyram2022direct  whose experimental data are published. Fig. 2b of that
work is an experiment performed using 200 fs pulses at 1040 nm, propagating
over 50 cm of GRIN fiber with 50 $\mu m$ diameter; 14 modal groups were
measured with differential eigenvalues $\epsilon_{i}$ ranging between 0 and
70671 m-1; thermalization was obtained in the experiment with 52 kW peak
power. The experimental mean modal power fractions $\lvert f_{i}\rvert^{2}$
are reported in Fig.7; degenerate modes are supposed with same modal power
fraction and eigenvalue.

Numerical fits were performed using the weighted BE, Eq. 1, and the RJ, Eq. 20
with $\gamma=1$. The RJ fit provided $T=400$ m-1 and $\mu^{\prime}=-72600$
m-1, equal to the values calculated in the original work; the RJ law appears
following the distribution up to the 7-th modal group, but it does not follow
correctly the distribution of the HOMs; the fit accuracy was $R^{2}=0.959$.
The weighted BE fit obtained $T=23300$ m-1, $\mu^{\prime}=-71100$ m-1 and
$\gamma=70.1$; the curve appears following the distribution up to the 14-th
group order, with accuracy $R^{2}=0.998$.

Figure 7: Experimental mean modal power fractions from Fig. 2b of
pourbeyram2022direct . RJ and weighted BE fits are reported.

Fig. 3e in pourbeyram2022direct  is a second experiment performed using 200 fs
pulses at 1040 nm, peak power 20 kW, propagating in 2.5 m of GRIN fiber with
62.5 $\mu m$ diameter. The measured modal groups were 28 with $\epsilon_{i}$
ranging between 0 and 176321 m-1; the measured mean modal power fractions are
reported in Fig.7. The RJ fit provided in this case $T=292$ m-1 and
$\mu^{\prime}=-178000$ m-1 with accuracy $R^{2}=0.938$; in the original work,
values were 300 and -178000 m-1 respectively; the RJ is able to follow the
experimental distribution up to the 8-th group order. The weighted BE fit
obtained $T=45700$ m-1, $\mu^{\prime}=-179000$ m-1 and $\gamma=166$; the fit
provided an accuracy of $R^{2}=0.976$ and it was able to overlap to the
experimental distribution up to the 17-th group order.

A validation of the state equation Eq. 2 was performed calculating the error
$\epsilon_{SE}$ from Eq. 3; the results were 0.64 % and 1.4 % for the two
experiments, respectively.

Figure 8: Experimental mean modal power fractions from Fig. 3e of
pourbeyram2022direct . RJ and weighted BE fits are reported.

## References

  *   * [1] Gloge, D. & Marcatili, E. A. J.  Multimode theory of graded-core fibers.  _The Bell System Technical Journal_ 52 (9), 1563–1578 (1973).  10.1002/j.1538-7305.1973.tb02033.x .
  * [2] Marcuse, D.  Derivation of coupled power equations.  _Bell System Technical Journal_ 51 (1), 229–237 (1972).  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/j.1538-7305.1972.tb01911.x.  https://doi.org/10.1002/j.1538-7305.1972.tb01911.x, https://onlinelibrary.wiley.com/doi/pdf/10.1002/j.1538-7305.1972.tb01911.x .
  * [3] Hasegawa, A.  Self-confinement of multimode optical pulse in a glass fiber.  _Opt. Lett._ 5 (10), 416–417 (1980).  URL https://opg.optica.org/ol/abstract.cfm?URI=ol-5-10-416.  10.1364/OL.5.000416 .
  * [4] Richardson, D. J., Fini, J. M. & Nelson, L. E.  Space-division multiplexing in optical fibres.  _Nature Photonics_ 7, 354 (2013).  10.1038/nphoton.2013.94 .
  * [5] Winzer, P. J., Neilson, D. T. & Chraplyvy, A. R.  Fiber-optic transmission and networking: the previous 20 and the next 20 years.  _Opt. Express_ 26 (18), 24190–24239 (2018).  URL https://opg.optica.org/oe/abstract.cfm?URI=oe-26-18-24190.  10.1364/OE.26.024190 .
  * [6] Wright, L. G., Christodoulides, D. N. & Wise, F. W.  Spatiotemporal mode-locking in multimode fiber lasers.  _Science_ 358 (6359), 94–97 (2017) .
  * [7] Wu, F. O., Hassan, A. U. & Christodoulides, D. N.  Thermodynamic theory of highly multimoded nonlinear optical systems.  _Nat. Photonics_ 13 (11), 776–782 (2019) .
  * [8] Fusaro, A., Garnier, J., Krupa, K., Millot, G. & Picozzi, A.  Dramatic acceleration of wave condensation mediated by disorder in multimode fibers.  _Phys. Rev. Lett._ 122, 123902 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.123902.  10.1103/PhysRevLett.122.123902 .
  * [9] Gloge, D.  Optical power flow in multimode fibers.  _The Bell System Technical Journal_ 51 (8), 1767–1783 (1972).  10.1002/j.1538-7305.1972.tb02682.x .
  * [10] Marcuse, D.  Losses and impulse response of a parabolic index fiber with random bends.  _Bell System Technical Journal_ 52, 1423–1437 (1973) .
  * [11] Ho, K.-P. & Kahn, J. M.  Linear propagation effects in mode-division multiplexing systems.  _J. Lightwave Technol._ 32 (4), 614–628 (2014).  URL https://opg.optica.org/jlt/abstract.cfm?URI=jlt-32-4-614 .
  * [12] Savović, S. _et al._ Power flow in graded-index plastic optical fibers.  _Journal of Lightwave Technology_ 37, 4985–4990 (2019) .
  * [13] Olshansky, R.  Mode coupling effects in graded-index optical fibers.  _Appl. Opt._ 14 (4), 935–945 (1975).  URL https://opg.optica.org/ao/abstract.cfm?URI=ao-14-4-935.  10.1364/AO.14.000935 .
  * [14] Zitelli, M., Mangini, F., Ferraro, M., Sidelnikov, O. & Wabnitz, S.  Conditions for walk-off soliton generation in a multimode fiber.  _Communications Physics_ 4 (1), 1–6 (2021).  URL http://dx.doi.org/10.1038/s42005-021-00687-0.  10.1038/s42005-021-00687-0 .
  * [15] Zitelli, M., Ferraro, M., Mangini, F. & Wabnitz, S.  Single-mode spatiotemporal soliton attractor in multimode GRIN fibers.  _Photonics Research_ 9 (5), 741 (2021).  10.1364/prj.419235 .
  * [16] Gordon, J. P.  Theory of the soliton self-frequency shift.  _Optics Letters_ 11 (10), 662 (1986).  10.1364/ol.11.000662 .
  * [17] Zitelli, M., Ferraro, M., Mangini, F. & Wabnitz, S.  Characterization of multimode soliton self-frequency shift.  _Journal of Lightwave Technology_ 40 (24), 7914–7921 (2022).  10.1109/JLT.2022.3205846 .
  * [18] Conti, C., Leonetti, M., Fratalocchi, A., Angelani, L. & Ruocco, G.  Condensation in disordered lasers: Theory, $3\mathrm{D}+1$ simulations, and experiments.  _Phys. Rev. Lett._ 101, 143901 (2008).  URL https://link.aps.org/doi/10.1103/PhysRevLett.101.143901.  10.1103/PhysRevLett.101.143901 .
  * [19] Klaers, J., Schmitt, J., Damm, T., Vewinger, F. & Weitz, M.  Bose–einstein condensation of paraxial light.  _Applied Physics B_ 105 (17) (2011).  URL https://doi.org/10.1007/s00340-011-4734-6.  10.1007/s00340-011-4734-6 .
  * [20] Parisi, G., Urbani, P. & Zamponi, F.  _Thermodynamics of Glass States_ , 99–139 (Cambridge University Press, 2020).
  * [21] Conti, C. & Leuzzi, L.  Complexity of waves in nonlinear disordered media.  _Phys. Rev. B_ 83, 134204 (2011).  URL https://link.aps.org/doi/10.1103/PhysRevB.83.134204.  10.1103/PhysRevB.83.134204 .
  * [22] Pourbeyram, H. _et al._ Direct observations of thermalization to a rayleigh–jeans distribution in multimode optical fibres.  _Nature Physics_ 18, 685–690 (2022) .
  * [23] Zitelli, M. _et al._ Spatiotemporal mode decomposition of ultrashort pulses in linear and nonlinear graded-index multimode fibers.  _Photon. Res._ 11 (5), 750–756 (2023).  URL https://opg.optica.org/prj/abstract.cfm?URI=prj-11-5-750.  10.1364/PRJ.484271 .
  * [24] Renninger, W. H. & Wise, F. W.  Optical solitons in graded-index multimode fibres.  _Nature Communications_ 4, 1–6 (2013).  10.1038/ncomms2739, arXiv:1302.1820 .
  * [25] Wright, L. G., Renninger, W. H., Christodoulides, D. N. & Wise, F. W.  Spatiotemporal dynamics of multimode optical solitons.  _Opt. Express_ 23 (3), 3492–3506 (2015).  URL https://opg.optica.org/oe/abstract.cfm?URI=oe-23-3-3492.  10.1364/OE.23.003492 .
  * [26] Krupa, K. _et al._ Spatial beam self-cleaning in multimode fibres.  _Nature Photonics_ 11 (4), 237–241 (2017) .
  * [27] Poletti, F. & Horak, P.  Description of ultrashort pulse propagation in multimode optical fibers.  _Journal of the Optical Society of America B_ 25 (10), 1645 (2008).  10.1364/josab.25.001645 .
  * [28] Berti, N. _et al._ Interplay of thermalization and strong disorder: Wave turbulence theory, numerical simulations, and experiments in multimode optical fibers.  _Phys. Rev. Lett._ 129, 063901 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevLett.129.063901.  10.1103/PhysRevLett.129.063901 .
  * [29] Jung, P. S., Wu, F. O., Hassan, A. U. & Christodoulides, D. N.  Measuring the optical temperature of a soliton.  _2019 Conference on Lasers and Electro-Optics (CLEO)_ 1–2 (2019) .
  * [30] Wu, F. O., Hassan, A. U. & Christodoulides, D. N.  Thermodynamic theory of highly multimoded nonlinear optical systems.  _Nature Photonics_ 13 (11), 776–782 (2019).  URL http://dx.doi.org/10.1038/s41566-019-0501-8.  10.1038/s41566-019-0501-8 .
  * [31] Podivilov, E. _et al._ Thermalization of orbital angular momentum beams in multimode optical fibers.  _Phys. Rev. Lett._ 128 (24), 243901 (2022) .
  * [32] Mangini, F. _et al._ Statistical mechanics of beam self-cleaning in grin multimode optical fibers.  _Optics Express_ 30 (7), 10850–10865 (2022) .
  * [33] Karlsson, M., Anderson, D. & Desaix, M.  Dynamics of self-focusing and self-phase modulation in a parabolic index optical fiber.  _Opt. Lett._ 17 (1), 22–24 (1992).  URL https://opg.optica.org/ol/abstract.cfm?URI=ol-17-1-22.  10.1364/OL.17.000022 .
  * [34] Ahsan, A. S. & Agrawal, G. P.  Spatio-temporal enhancement of raman-induced frequency shifts in graded-index multimode fibers.  _Opt. Lett._ 44 (11), 2637–2640 (2019).  URL https://opg.optica.org/ol/abstract.cfm?URI=ol-44-11-2637.  10.1364/OL.44.002637 .
  * [35] Ahsan, A. S. & Agrawal, G. P.  Graded-index solitons in multimode fibers.  _Opt. Lett._ 43 (14), 3345–3348 (2018).  URL https://opg.optica.org/ol/abstract.cfm?URI=ol-43-14-3345.  10.1364/OL.43.003345 .
  * [36] Gloge, D.  Bending loss in multimode fibers with graded and ungraded core index.  _Appl. Opt._ 11 (11), 2506–2513 (1972).  URL https://opg.optica.org/ao/abstract.cfm?URI=ao-11-11-2506.  10.1364/AO.11.002506 .
  * [37] Flamm, D., Naidoo, D., Schulze, C., Forbes, A. & Duparré, M.  Mode analysis with a spatial light modulator as a correlation filter.  _Optics Letters_ 37 (13), 2478 (2012).  10.1364/ol.37.002478 .
  * [38] Jaynes, E. T.  Gibbs vs boltzmann entropies.  _American Journal of Physics_ 33 (5), 391–398 (1965).  URL https://doi.org/10.1119/1.1971557.  10.1119/1.1971557, https://doi.org/10.1119/1.1971557 .
  * [39] Poletti, F. & Horak, P.  Description of ultrashort pulse propagation in multimode optical fibers.  _J. Opt. Soc. Am. B_ 25 (10), 1645–1654 (2008).  URL https://opg.optica.org/josab/abstract.cfm?URI=josab-25-10-1645.  10.1364/JOSAB.25.001645 .
  * [40] Wright, L. G. _et al._ Multimode nonlinear fiber optics: Massively parallel numerical solver, tutorial, and outlook.  _IEEE Journal of Selected Topics in Quantum Electronics_ 24 (3), 1–16 (2018).  10.1109/JSTQE.2017.2779749, arXiv:1708.05324 .