Statistics of modal condensation in nonlinear multimode fibers

Optical pulses traveling through multimode optical fibers encounter the influence of both linear disturbances and nonlinearity, resulting in a complex and chaotic redistribution of power among different modes. In our research, we explore the phenomenon where multimode fibers reach stable states marked by the concentration of energy into both single and multiple sub-systems. We introduce a weighted Bose-Einstein law, demonstrating its suitability in describing thermalized modal power distributions in the nonlinear regime, as well as steady-state distributions in the linear regime. We apply the law to experimental results and numerical simulations. Our findings reveal that, at power levels situated between the linear and soliton regimes, energy concentration occurs locally within higher-order modal groups before transitioning to global concentration in the fundamental mode within the soliton regime. This research broadens the application of thermodynamic principles to multimode fibers, uncovering previously unexplored optical states that exhibit characteristics akin to optical glass.


Introduction
Multimode (MM) optical fibers [1][2][3] 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) [4,5], and the possibility of up-scaling the pulse energy delivered from fiber lasers [6].
In recent years, a thermodynamic interpretation of beam propagation in multimode optical fibers and systems has been formulated [7,8].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 β 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) [9][10][11], 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 [9,12,13], 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 FWHM depending on the pulse wavelength [14], and with energy given by E s = 1.76λ|β 2 (λ)|w 2  e /(n 2 T FWHM ), with β 2 (λ) (s 2 /m) the chromatic dispersion, n 2 m 2 /W the Kerr nonlinear coefficient, and w e the modal effective waist.At telecom wavelength λ = 1550 nm, the soliton pulsewidth is T 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 [15], 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 [16,17] (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 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 [18] and Bose-Einstein condensation [19].
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 [20,21]; this condition appears as intermediate between the low-energy disordered state, and the high-energy, highly condensed soliton state.

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 10 5 to 10 9 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) In Eq. 1, |f i | 2 = 2n j /(γn 0 g i ) is the mean modal power fraction over two polarizations, with γ given by the total number of energy packet N = γ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], i = β i − β j=Q are the differential modal eigenvalues with β i (m −1 ) the propagation constants.Eq. 1 was obtained under the reasonable condition n 0 >> exp − µ + i T − 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 µ (m −1 ), temperature T (m −1 ) and the normalized power γ are three degrees of freedom for fitting Eq. 1 to the experimental data, with the only constraint for γ to scale with the input power, and to respect the conservation law 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 an experimental error upon the state equation can be calculated as In Eq. 2, the first term in the sum (the normalized energy) must be constant; µ and T /γ must scale in order to maintain constant the remaining part of the equation.Not necessarily the parameters T /γ and µ 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 |f Since the input modal distribution is unchanged at high power, Eq. 2 can be used, together with |f (out) j | 2 measured at low power, to predict the thermodynamic parameters of the system at high power, when thermalization, or condensation, is achieved [22].However, in long spans of fiber at low power, RMC causes |f | 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) [16] 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 (β 1 − β 10 )/β 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.

Experiments
In order to determine the mode power distribution at the output of GRIN MM fibers, we used the mode-decomposition method introduced in [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 µm and a 10 µ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 nearfield (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 ps 2 /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).Quasisoliton 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 dispersioninduced temporal broadening) is as short as L ef f = 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 [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 |f i | 2 = 2P j /(g i P tot ), by supposing power equipartition within each group (see Methods Sec 6.3).
Fig. 2a shows that, in the linear propagation regime (i.e., for E in ≤ 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 [14]; 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 ≥ 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 [16,17,24,25].
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 [26].
Thanks to the accuracy offered by the used decomposition method, Fig. 3 reports in log scale the average output power fractions |f i | 2 from Fig. 2, versus their respective mode eigenvalues 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 [26].Output modal distributions are fitted using Eq. 1, taking care to respect the constraint for the normalized power γ.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 [12] [9] (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 µ = −58.2mm −1 , T = 8.71 mm −1 and γ = 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 ≥ 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 µ = −53.1(−52.Table 1 shows the values of the optical temperature T , differential potential µ , factor γ and state equation error 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, µ reduces from -65.42 mm −1 to -52.08 mm −1 ; γ 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.A detailed information on how the relative power fraction for groups of nondegenerate modes, or j • |f j | 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 [22].

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 [12] [9] (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 α 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 × 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 |f i | 2 vs. the mode eigenvalue progressively increases (Fig. 5b).This distribution can be fitted well by the weighted BE law with µ = −54.6 mm −1 , T = 31.1 mm −1 and γ = 70.5, which describes the RMC-induced steady-state.

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 coupledmode GNLSEs [27], 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 [12], 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.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 • |f j | 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 [15].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.

Discussion and Conclusions
The mean modal power fractions |f i | 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.Powerflow 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, quasisoliton 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 [28]; 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 [14].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 [17].
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 [29].All of the observed distributions correspond to steady states, being their modal distributions properly fitted by the weighted BE.

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 [30], 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 System entropy is defined as S = ln(W ); by applying the Stirling approximation, valid for n j + g j − 1 >> 1 we obtain (n j +g j −1) ln(n j +g j −1)−1 −(g j −1) ln(g j −1)−1 −n j ln(n j )−1 ; (5) Global system thermalization corresponds to the extremization of Eq. 5 using Lagrange multipliers for the total number of particles N = j n j and total normalized energy E = j β j n j , being β j the modal propagation constants supposed equal for the degenerate modes which provides, with no approximations Eq. 7 is valid for We can choose a = µ/(T n 0 ) and b = 1/(T n 0 ), with n 0 a reference number of energy packets (for example the value at the lowest tested power) and N = γn 0 .T (1/m) is an optical temperature and µ (1/m) a chemical potential.We these choices, Eq. 8 can be written as An alternative development of Eq. 8 consists in replacing a and b with non-factorizable constants a and b defined as The approximation in Eq. 10 is valid for |a + bβ j | << 1, which is less than 10 −5 in the experiments.By replacing into Eq.8 we obtain The particular choice of a, b, a and b conserves the system's energy and power; in fact it results, for a , b and For a, b it also results hence, the choices of a, b or a , b are equivalent in terms of power and energy conservation.
When passing from Eq. 9 to Eq. 11 by replacing the a , b 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 µ and T are now related to the particular set of β 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.
Introducing the differential eigenvalues i = β i − β j=Q , referred to the higher mode, and µ = µ + β j=Q , eq. 11 provides the weighted Bose-Einstein (BE) modal distribution The constraints for γ are to scale with the input power, and to respect the conservation law 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, γ = N/n 0 = P/P 0 .The system's internal energy is U = − j β j n j P/N (W/m); the power P = j n j P/N (W).From Eqs. 14 and 15 we get which provides the state equation From the local extremization problem, providing Eq. 11, it is not possible to derive the well-known Rayleigh-Jeans (RJ) distribution [30] [31] under reasonable approximations.However, from the global extremization solution, Eq. 8 with the definitions of a and b , we obtain the presence of N = γ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 |µ + i | << |T n0|, which brings directly to the RJ distribution In terms of fractional power , by choosing T = P 0 T (W/m), Eq. 20 provides.
The RJ law demonstrates to be a good choice for describing experiments characterized by global condensation states, such as the well-known selfcleaning experiments [26,31,32], while the more general weighted BE, Eq. 16, is suitable to fit the experimental data also in cases where local condensed states are achieved.

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 λ/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 µm.The corresponding selfimaging induced beam compression factor C = 2z p /(πβ 0 w 2 0 ) = 0.305, where β 0 = 2πn 0 /λ, z p = 0.55 mm the self-imaging period, and n 0 = 1.46 is the core index [33,34]; the effective beam waist for determining the nonlinear coefficient is w e = √ Cw 0 = 7.27µm [34,35], 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 µm, in order to increase the contribution from higher-order modes.
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 ∆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 √ z) [10].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 µW resolution.
Traditional 2D modal decomposition methods [37] 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 [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).

Power-flow numerical model
In the linear regime, RMC is properly modelled by the well-known powerflow diffusive equations.If P j is the power of the j-th mode group, the power exchange among adjacent modal groups is described by [12] with m(j) = j − 1, D and α 0 (1/m) the coupling coefficient and linear losses, respectively, A the modal loss coefficient, ∆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 |f i | 2 = 2P j /(g i P tot ).

Modal power equipartition
It is possible to demonstrate, using the alternative Gibb's definition of entropy [38], that modal equipartition does not generally apply among different groups.Given, at steady state, p i =< |f i | 2 > the modal occupation probability, and λ 1 , λ 2 two Lagrange multipliers, extremization of the entropy, while power P and energy E are conserved, reads p i β i − E/P = 0; (23) from Eq. 23, if only the power is conserved (λ 2 = 0), we obtain which brings to the modal equipartition p i = 1/M .However, when considering also the conservation of the energy, Eq. 23 provides which provides equipartition only for degenerate modes, with same β i .

GNLSE numerical model
In the nonlinear (and linear) regime, numerical simulations use the coupled GNLSEs [39], modified to include modal and wavelength-dependent losses, and linear random-mode coupling (RMC) for mode p In Eq. 26, β (p) n is the n-th order dispersion term (modal and chromatic) for mode p, α p (λ) the modal and wavelength-dependent loss coefficient, n 2 (m 2 /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 [12] q mp = 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 λ = 1400 nm are: β 2 = −11.8ps 2 /km, β 3 = 0.102 ps 3 /km for the fundamental mode, α = 2.6 dB/km, negligible modal loss, n 2 = 2.7x10 −20 m 2 /W, D = 0.003 (1/m), L c = 6 mm, f R = 0.18, S plmn calculated from modal overlap integrals. ------- The RJ fit provided in this case T = 292 m −1 and µ = −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 8th group order.The weighted BE fit obtained T = 45700 m −1 , µ = −179000 m −1 and γ = 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 SE from Eq. 3; the results were 0.64 % and 1.4 % for the two experiments, respectively.
Fig. A2 Experimental mean modal power fractions from Fig. 3e of [22].RJ and weighted BE fits are reported.

Fig. 1
Fig. 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.

Fig. 2
Fig. 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. 3
Fig. 3 Average mode power fraction |f i | 2 vs. modal eigenvalue 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.
2) mm −1 , T = 8.71(10.14)mm −1 and γ = 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 (µ + 1 )/T −0.13, which means that the RJ approximation to the BE law (i.e., |f i | 2 ∝ −T /(µ + i )) is appropriate around the fundamental mode with i = 1.

Fig. 4 a
Fig. 4 a) Experimental group power fraction j • |f j | 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.
Figure also shows that a simple exponential function, of the type |f i | 2 = a exp (−b i ) cannot properly fit the HOM groups for low losses.

Fig. 5 a
Fig. 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.

Fig. 6
Fig. 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
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;

q
mp A m + in 2 k 0 l,m,n

1 DLc
from modes m with g m = g p − from modes m with g m = g p + 1,