Open Access
Add to BibSonomy
Abstract
Since its first demonstration in graded-index multimode fibers, spatial beam self-cleaning has attracted a growing research interest. It allows for the propagation of beams with a bell-shaped spatial profile, thus enabling the use of multimode fibers for several applications, from biomedical imaging to high-power beam delivery. So far, beam self-cleaning has been experimentally studied under several different experimental conditions. Whereas it has been theoretically described as the irreversible energy transfer from high-order modes towards the fundamental mode, in analogy with a beam condensation mechanism. Here, we provide a comprehensive theoretical description of beam self-cleaning, by means of a semi-classical statistical mechanics model of wave thermalization. This approach is confirmed by an extensive experimental characterization, based on a holographic mode decomposition technique, employing laser pulses with temporal durations ranging from femtoseconds up to nanoseconds. An excellent agreement between theory and experiments is found, which demonstrates that beam self-cleaning can be fully described in terms of the basic conservation laws of statistical mechanics.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Interest in optical multimode fibers (MMFs) has experienced different periods of rise and fall, mostly driven by their pros and cons in several applications, particularly in optical communication systems [1]. On the one side, MMFs are capable of propagating multiple transverse modes, which enables a substantial increase of the information capacity of fiber-optic links, based on the mode-division-multiplexing technique [2]. Furthermore, thanks to their large core area, MMFs deliver higher energies with respect to singlemode fibers, which is advantageous for developing novel high-power fiber laser architectures [3]. On the other hand, large modal temporal dispersion has prevented the use of MMFs in long-distance high-bit-rate transmissions. Moreover, since the multitude of initially excited modes spatially overlaps, the beam quality at the output of MMFs is severely degraded by multimode interference (or modal noise), when compared to singlemode fibers. As a result, nowadays singlemode fibers are almost exclusively employed.
In recent years, interest in the use of MMFs has returned within coherent optical communication systems, thanks to their capability for pre- or post-compensation of linear modal dispersion effects [2,4,5]. Moreover, several novel nonlinear optical effects have been explored [6–9]. Remarkably, it has been shown that the Kerr effect can be usefully exploited, in order to provide a bell-shape to the beam profile at the output of MMFs, which permits a significant beam quality and brightness improvement [10]. This effect has been dubbed as spatial beam self-cleaning (BSC) [11]. BSC has been demonstrated in several configurations, either in the normal [10–21] or in the anomalous dispersion regime [22–24], in graded-index (GRIN) as well as in step-index fibers [25], and with input pulse durations ranging from nanoseconds to femtoseconds.
Many of these studies have investigated the physical mechanism behind BSC. It is generally accepted that BSC relies on nonlinear coupling among the multiple modes of GRIN MMFs, which leads to a modal redistribution of the input beam energy [26]. Moreover, the robustness of the BSC effect indicates that a sort of equilibrium distribution is established as a result of beam propagation in the MMF. This has led to a description of BSC in a thermodynamics framework. The bell-shaped beam that is spontaneously generated suggests that a prevailing population of the fundamental mode is reached. However, the size of self-cleaned beams was measured to be generally wider than that of the fundamental mode of the MMF [8]. In addition, the output bell-shaped beam is typically sitting on a wide low-power background. This indicates that higher-order modes (HOMs) also contribute to determining the observed output beam shapes. In this context, it was pointed out that BSC can be seen as a thermalization phenomenon, where the pair ($T, \mu$), i.e., temperature and chemical potential, is fixed by the pair ($N,H$), i.e., number of particles and internal energy [27–29].
Given that BSC manifests itself when the power of the laser beam that is injected into the fiber overcomes a certain threshold value, one might think of BSC as a sort of (quantum) Bose-Einstein condensation of classical waves [30]. This is only partially true: as a matter of fact, BSC can be associated to a wave condensation phenomenon, as Baudin and coworkers have recently demonstrated in [31]. However, in that work, it has been shown that at a given input laser power, BSC can be achieved when the temperature $T$ drops below a certain critical value, i.e., only when the input coupling condition (CC), which determines the ratio $h=H/N$, is favorable enough (cfr, the horizontal arrow in Fig. 1). Such a phase transition arises when $h$ decreases below a critical value, say, $h_c$ [30].
Fig. 1. Optical phase diagram illustrating orthogonal paths to describe BSC as a thermalization or a condensation process, respectively. The diagram is split into two zones. Below the solid black line (which indicates the condition for observing BSC), optical beams are in a disordered phase (speckled beams). Whereas self-cleaned beams only exist in the upper-right region of energy per particle $h$ (or coupling conditions, $CC$) and number of particles (or input power $P$). When passing from one zone to the other by a vertical path (red arrow), i.e., by varying $N$ while keeping constant $h$, BSC can be described as a thermalization of fiber modes. The latter follows the RJ distribution, and HOMs have a non-negligible population (white bars indicate the equilibrium modal occupancy). On the other hand, when keeping fixed $N=N_0$ and decreasing $h$ below its critical value (horizontal phase transition or grey arrow), BSC can be described as a wave condensation phenomenon, i.e., only the fundamental mode has a macroscopic population.
On the other hand, in the thermodynamic description the input laser power plays the role of the number of particles ($N$). Therefore, increasing the input laser power does not necessarily trigger wave condensation. As a matter of fact, a thermodynamic system can reach thermal equilibrium even at temperatures which are higher than the condensation temperature, i.e., with CC corresponding to $h>h_c$. Of course, the input power must be sufficiently high for allowing the occurrence of significant degenerate four-wave-mixing (FWM) among the modes. In other words, $N$ must be large enough ($N>N_c$, solid black line in Fig. 1) to reach thermal equilibrium over a finite propagation distance (ergodicity condition [27]).
Whenever BSC is observed by increasing the input power (or $N$) for a given input CC (independently on whether $h$ is greater or smaller than $h_c$), the effect can be described as an example of wave thermalization (cfr, the vertical arrow in Fig. 1): here we follow this approach. At thermodynamic equilibrium, the fiber mode distribution obeys the Rayleigh-Jeans (RJ) law, in agreement with the thermodynamic analysis of classical systems comprising a finite number of modes [27]. In a nutshell, thermalization of the mode distribution is a necessary condition for condensation, but not viceversa. As sketched in Fig. 1, the mode distributions corresponding to wave condensation or thermalization are generally different. Because the former involves a dominant (or macroscopic) population of the fundamental mode only, accompanied by energy equipartition among all HOMs. Whereas the latter is characterized by a RJ mode distribution that, although the fundamental mode is always the most populated, may also contain a macroscopic contribution from the adjacent low-order modes.
Note that, by properly adjusting the input beam CC or by means of wavefront shaping, BSC into the $LP_{11}$ [32] or even several HOMs [33] could be observed. However, experiments and numerical simulations reported in the supplementary of [21] indicate that such HOMs are transient states, which ultimately decay into a mode distribution dominated by the fundamental mode. In these cases, the input beam power that is required for thermalization can be higher than the threshold of competing nonlinear effects, such as intermodal FWM and stimulated Raman scattering, which deplete the input beam.
From an experimental point of view, a thermodynamic analysis of BSC requires the use of mode decomposition (MD) techniques for a proper analysis of the output beam mode populations. Several MD methods have been proposed in the literature. Conceptually, they can be divided into different groups, depending on the foundation lying beneath them: genetic algorithms [34], adaptive optics [35], digital holography with spatial light modulators (SLM) [36], mode-cross correlation analysis [37], and beam characterization in terms of spatiotemporal [38] or spatial-spectral [39] distributions.
In this work, we present a comprehensive and yet simple description of the phenomenon of BSC within the thermodynamic framework. We derive a mathematical model starting from statistical mechanics, pointing out the dependence of the equilibrium RJ distribution on the radial and azimuthal properties of GRIN fiber modes. In particular, we point out that, if the modes are properly sorted, the RJ distribution can be described as a sequence of sub-equilibrium distributions. It turns out that the modes which have radial symmetry are the most populated in the occurrence of BSC. Whereas, the mode occupancy vanishes when the spatial (azimuthal) complexity grows. Theoretical predictions are compared to an extensive MD experimental study, based on the holographic method [40]. We investigate BSC when operating with pulses ranging from the femto- to the pico- and nanosecond regime. Our results show that BSC can be reached independently of the input pulse duration, and that the parameters of the output RJ distribution can be tuned by acting on the laser-fiber injection condition, or CC. The obtained reconstructions of the spatial profile of the beam at the fiber output turn out to be remarkably similar to the measured patterns, which proves the power and the accuracy of our MD method. Finally, we experimentally verify the statistical mechanics conservation laws, which are at the basis of our theoretical model.
2. Theory
Here, we propose a semi-classical approach to determine the mode equilibrium distribution at the occurrence of BSC. In GRIN fibers, modes are conventionally defined in terms of a pair of integer numbers ($\ell,m$), which indicate their radial and azimuthal properties, respectively. As a remarkable example, we may consider Laguerre-Gauss functions as the modal basis, which will be the case for our MD experiments (see Appendix Section). In our derivation, we also refer to the quantum number $q = 2\ell + |m|$, so that $m = -q, -q+2,\ldots, q-2,q$. An illustration of GRIN fiber modes and their relative sorting is provided in Fig. 2(a). In Fig. 2(b), we report the couple ($\ell$,m) which is associated to each mode in Fig. 2(a). The upper limit of $q$, which we dub $Q$, is defined by the fiber cut-off, i.e. by the MMF geometrical parameters and the laser wavelength [41]. The presence of a frequency cut-off is a fundamental condition for avoiding the divergence of the entropy at low temperatures, which is referred to as ultraviolet catastrophe in gas thermodynamics [30]. For the sake of simplicity, in Fig. 2(a) we show modes up to $q = 5$. The dynamics of the system of modes can be described by means of a Hamiltonian operator, which only takes into account the linear kinetic term:
where $p_{\ell,m}$ and $n_{\ell,m}$ are the momentum and the occupancy of mode ($\ell,m$), and $c$ is the speed of light in vacuum. Since it has been shown that BSC mainly involves a purely spatial dynamics [11,19], we neglect all temporal effects, and consider monochromatic waves only.Fig. 2. a) Illustration of the GRIN MMF modes displayed in the ($q,m$) plane. b) Relative mapping of the indexes ($\ell,m$). c) GRIN fiber modes momentum vs. modes index sorted by quantum number $q = 2\ell +|m|$. d) RJ distribution at the equilibrium in the ($q,m$) plane, when setting $N = 1$, $\mu =-60$ $\mathrm{mm}^{-1}$, and $T=0.3$ $\mathrm{mm}^{-1}$. The values of $k_{0,0}$ and $k_\text {SSI}$ are calculated starting from the nominal values of numerical aperture and relative core-cladding index difference provided by the fiber manufacturer.
In the framework of statistical mechanics, we impose that the total energy of the system and the number of particles are conserved upon propagation, i.e., the Hamiltonian (1) is constant and
The condition (2) holds as long as any dissipative effects can be neglected, which is always the case in practical demonstrations of BSC in passive MMFs. We underline that the nonlinear dynamics of mode coupling in MMFs should lead to a proper description in the framework of the micro-canonical ensemble, i.e., the statistical ensemble in which both $N$ and $T$ are conserved. However, the derivation of the equilibrium distribution in this ensemble may be rather challenging. For this reason, here we rely on the Gibbs theorem on the equivalence of ensembles, which states that, in the thermodynamic limit, the equilibrium distribution must the same in all of the statistical ensembles [42]. Therefore, for the sake of simplicity, we develop our theory in the frame of the grand canonical ensemble, where the temperature and the chemical potential are fixed, at variance with the energy and the number of particles. Specifically, dubbing $\beta$ and $\mu$ the Lagrange’s multipliers corresponding to each conservation law, the partition function reads [43]:
Equation (4) is the RJ distribution, which is reached at thermal equilibrium under the hypothesis of the ergodic theorem. The latter holds whenever $N$ is sufficiently large for ensuring a complete exchange of energy among the modes via FWM processes, so that the multimode system explores its whole phase space. Experimentally, this condition reflects the presence of a threshold of the input laser power for achieving BSC. It is important to underline that the ergodicity condition can also be matched, at least in principle, by extending the evolution ”time” of the system, i.e., by increasing the fiber length. It is well-known, in fact, that the power threshold for BSC decreases as the fiber length increases [11]. Very often, in fact, people tune the input power as an equivalent cut-back experiment when dealing with nonlinear mode coupling, e.g., when studying soliton self-mode conversion in MMFs [44]. Nonlinear mode coupling is particularly effective in GRIN fibers, which, thanks to their parabolic profile of the core refractive index, force the beam to periodically shrink after a propagation distance $z_\text {SSI}$, where SSI stands for spatial-self imaging [45]. This peculiarity originates from the fact that the modes of GRIN fibers have equispaced momenta, i.e.,
where $k_{0,0}$ is the propagation constant of the fundamental mode, which can be calculated by starting from the numerical aperture and the core size (usually provided by the fiber manufacturer), andIn Fig. 2(c), we display the modes’ momenta, where the modes are sorted by the quantum number $q$ by following the lexicographical order with respect to Fig. 2(b). Finally, one can calculate the average number of particles in mode ($\ell,m$) as
where we defined $\beta = 1/K_B T$, as it is conventionally done in thermodynamic systems, where $K_B$ is the Boltzmann’s constant. Moreover, $\hbar$ is the Planck’s constant, while $T$ and $\mu$ are parameters which only depend on the input CC, and need to be determined by MD experiments. As $T$ and $\mu$ are only statistical quantities, in the following we will set $K_B = \hbar = c = 1$. In this way, $T$, $\mu$ and $H$ have the dimension of the inverse of a length. In Fig. 2(d), we plot the thermal distribution (7) in the ($q,m$) space. Note that a correspondence between the thermodynamical variables of classical gases and multimode gases has been recently reported in Ref. [46].3. Experimental setup
The MD experimental setup that we used to study BSC in GRIN MMFs is shown in Fig. 3. It consists of an ultra-short pulse laser system pumped by a femtosecond Yb-based laser (Light Conversion PHAROS-SP-HP), generating pulses with adjustable duration (by means of a dispersive pulse stretcher), at 100 kHz repetition rate and $\lambda =1030$ nm, and with Gaussian beam shape ($M^{2}$=1.3). The pulse shape was measured by using an autocorrelator (APE PulseCheck type 2), resulting in a sech temporal shape with pulse widths ranging from 174 fs up to 10 ps. As shown in Fig. 3, the laser beam was injected by a convergent lens ($L_0$) into the core of the GRIN fiber. The input diameter at $1/e^{2}$ of peak intensity was measured to be $30$ $\mu$m. We employed 3 m long standard 50/125 GRIN fibers (GIF50E from Thorlabs), whose core radius, core refractive index along the axis, relative core-cladding index difference, numerical aperture and fundamental mode radius at $\lambda$ = 1030 nm are $r_c=25$ $\mu$m, $n_0=1.472$, $\Delta =0.0103$, $NA = 0.2$ and $r_{0,0} = 6.33$ $\mu$m, respectively. The near-field profile at the fiber output is imaged onto an SLM (Hamamatsu LCOS- X15213) by means of two confocal lenses ($L_1$, with $f_1=2.75$ mm and $L_2$, with $f_2=400$ mm). Between those, we placed a bandpass filter (BPF, $1030 \pm 5$ nm), a half-wave plate ($\lambda /2$), and a polarizer (PBS). Our measurement system allows us to avoid the parasitic influence of nonlinear frequency conversions, e.g., provided by Raman scattering or geometric parametric instability [10], which is detrimental for our MD reconstruction algorithm. As a matter of fact, the phase pattern on the SLM for the profile reconstruction algorithm must be chosen at a given wavelength. We could also tune the intensity of the beam reaching the SLM by means of the $\lambda /2$ wave plate. It is worth to underline that selecting an output polarization does not invalidate our theoretical description of BSC. According to Ref. [47], if enough input power is provided, BSC occurs independently of the output polarization state. A flip mirror (FM) is used for imaging the near-field profile at the fiber output facet onto an IR camera (NF, Gentec Beamage-4M-IR). Images acquired in this way were used as a reference, in order to check the quality of the reconstruction made by the MD algorithm. At last, a convex lens ($L_3 = 400$ mm) projects the field reflected by the SLM onto a second camera (FF camera). The lens is placed in the middle between the SLM and the camera, so that both of these objects are at its focal distance. Finally, the beam average power at the input and the output of the fiber was measured by a photodiode power meter (Thorlabs). In order to explore the nanosecond pulse regime, we carried out analogous MD experiments by using a Nd:YAG laser emitting pulses of 0.435 ns duration and 1064 nm of wavelength, at 1 kHz of repetition rate.
Fig. 3. Sketch of the experimental setup. The lenses focal distances are $f_1 = 2.75$ mm, and $f_2 = f_3 = 400$ mm. The arrows indicate the beam polarization direction (horizontal).
4. Results
Since the values of the intensity profiles which are detected by the cameras are given in arbitrary units, our setup only allows for evaluating the relative values of mode occupancy, i.e., the power fraction values, normalized by the total beam power. Therefore, in the following of our analysis we set $N=1$, so that $n_{\ell,m}$ represents the occupancy fraction of mode ($\ell,m$). Moreover, we limited our experimental analysis to the first 78 modes with the highest values of momentum, i.e., we only considered modes with $q\leq 11$ (see Fig. 2(c)). We found out that such a number was a good compromise between the time convergence of the MD algorithm, and the quality of the near-field reconstructions. Further details about the MD method are presented in the Appendix section.
In Fig. 4, we report results obtained when employing 7.6 ps input laser pulses. Specifically, what we show in Fig. 4(a) is a 3D histogram of the experimental mode power fractions, for different values of the input beam peak power ($P_p$). Whereas, in Fig. 4(b), we illustrate the experimental evolution, as a function of $P_p$, of the sum of the power fraction of all modes with the same quantum number. This figure shows that the occupancy of the fundamental mode (which is the only mode with $q=0$) grows larger when $P_p$ is increased, indicating the progressive occurrence of BSC, accompanied by thermalization of the mode distribution. In particular, it can be seen that at the lowest input peak power value of 110 W, the CC is such that the sum of the occupancy of modes with $q=1$ is higher than that of the fundamental mode. The latter progressively becomes dominant, when increasing the input power, up to the occurrence of BSC. In Fig. 4(c), we compare the experimental power fraction values (blue bars) with the theoretical RJ distributions (red dashed lines). As it can be seen, as $P_p$ increases, the experimental distribution progressively approaches the theoretical RJ distribution. This can be quantitatively appreciated in Fig. 4(d), where we plot the root-mean-square error (RMSE) of the observed mode distributions with respect to the RJ distribution, as a function of $P_p$. Finally, insets of Fig. 4(c) demonstrate the excellent quality of our MD method. Here, we show that the measured output near-field intensities (images in the left column) are impressively similar to their MD reconstructions (images in the right column), for all input power values.
Fig. 4. MD analysis for 7.6 ps of input pulse duration. a) Mode distribution (with the same sorting as in Fig. 2(c)) for different values of the input peak power. b) Experimental value of the power fraction of the fiber modes grouped by $q$, vs. input power $P_p$. c) Details of the reconstruction in 2D plots. The images in the inset of each graph represent the measured (left) and reconstructed (right) output beam profile. The solid line is the RJ distribution that fits the experimental values at the highest value of input power (8.11 kW). The fitting parameters are $\mu = -63.97$ $\mathrm{mm}^{-1}$ and $T = 0.92$ $\mathrm{mm}^{-1}$. d) Root-mean-square error (RMSE) of the fitting curve in b), when varying $P_p$.
4.1 BSC stability vs input pulse duration
We carried out MD experiments with pulses whose durations range from hundreds of femtoseconds up to 0.5 nanoseconds. In all cases, we kept the same fiber length, whereas we used the Nd:Yag laser that emits pulses of 0.5 ns of duration. In Fig. 5, we present the MD analysis based on 174 fs input pulses. This figure allows for further appreciating the thermalization nature of BSC. We could observe, in fact, that the fundamental mode power fraction has a non-monotonic behavior when increasing $P_p$, although the RMSE with respect to the RJ distribution keeps decreasing. At $P_p \simeq 26$ kW, the power fraction $\langle n\rangle$ of the fundamental mode (with $q=0$) has a local minimum (see Fig. 5(a) and b). It is important to underline that MD experiments of Fig. 4 and Fig. 5 with different pulse duration were carried out with the same laser-fiber CC. Our laser allows, in fact, to vary the pulse duration, while the pulse energy is kept constant. Therefore, we could maintain the experimental setup untouched, thus avoiding any modification of the injection conditions.
Fig. 5. Same as Fig. 4 for pulses of 174 fs of duration. The CC is the same as in Fig. 4. The fitting parameters are $\mu = -63.96$ $\mathrm{mm}^{-1}$ and $T = 1.01$ $\mathrm{mm}^{-1}$.
We found almost identical values of the parameters $\mu \simeq 64$ $\mathrm{mm}^{-1}$ and $T\simeq 1$ $\mathrm{mm}^{-1}$ in these two cases. In agreement with theoretical predictions, in fact, $T$ and $\mu$ do not depend on the pulse duration, and they can be varied only by changing the CC [31]. This is the case of the MD analysis which is reported in Fig. 6 and 7, where we found different values of $T$ and $\mu$ with respect to those in Fig. 4 and 5. Besides their different laser-fiber CC, the results reported in Fig. 6 and 7 have been obtained with a different pulse duration of 1 picosecond, thus experimentally confirming that BSC can be achieved with pulses ranging from the femtosecond up to the sub-nanosecond regime.
Fig. 6. Same as Fig. 4 for pulses of 1 ps of duration, and a different CC from Figs. 4 and 5. The fit parameters are $\mu = -64.26$ $\mathrm{mm}^{-1}$ and $T = 1.08$ $\mathrm{mm}^{-1}$.
Fig. 7. Same as Fig. 4 for pulses of 435 ps of duration and different injection conditions. The fit parameters are $\mu = -51.86$ $\mathrm{mm}^{-1}$ and $T = 0.72$ $\mathrm{mm}^{-1}$.
4.2 Radial modes distribution
So far, when showing modal distributions, we sorted the modes by their quantum number. This is the most convenient choice when one limits the analysis to beams which do not carry angular momentum, i.e., when injecting a radially symmetric laser beam at the center of the fiber core, as it occurs in our case. However, our theoretical derivation aims to be more general, since it has been developed on the basis of Laguerre-Gauss modes, which depend on two indexes (radial and azimuthal). Then one may naturally wonder how the equilibrium distribution looks like when sorting the modes by $\ell$ and $m$. This is shown in Fig. 8, for all considered values of the input pulse duration.
Fig. 8. Illustration of the mode decomposition results when sorting the modes by their radial index $\ell$ for input pulse duration of 174 fs (a), 1 ps (b), 7.6 ps (c), and 435 ps (d), respectively. Black bars represent the experimental values of the power fraction of each mode, while the red dashed line is the theoretical RJ distribution.
Figure 8 permits to clearly appreciate that the RJ distribution is indeed a sequence of sub-distributions. As a matter of fact, for each value of the radial index $\ell$, the mode having $m=0$ is the most populated at the equilibrium. Whereas the occupancy $\langle n_{\ell,m\neq 0}\rangle$ rapidly vanishes as $m$ grows. We underline that the steps of the red solid curve (the theoretical RJ distribution) when passing from $\ell$ to $\ell +1$ are due to the finite number of modes of the system, which in our case is due to the mode truncation that we made when running our MD algorithm. However, similar curves, albeit much longer, would have been obtained, had we considered the presence of all modes below cut-off.
When modes are sorted by the radial index, their distribution emphasizes the symmetry with respect to the azimuthal index $m$ of the equilibrium distribution. From a mathematical point of view, this is provided by the presence of the absolute value in the definition of the momentum of the GRIN fiber modes, see Eq. (5). Thus, according to Eq. (7), $\langle n_{\ell,m}\rangle = \langle n_{\ell, -m}\rangle$. Experimentally, this is demonstrated by Fig. 8. As it can be seen, each couple of bars corresponding to modes with the same value of $\ell$, and opposite signs of $m$ have the same height, again independently of the pulse duration.
4.3 Conservation laws verification
Further investigating the symmetry of mode occupancies with respect to $m$ highlights an important feature of BSC: the conservation of the mode parity ($M$), which is defined as:
The conservation of $M$ is related to the conservation of the total angular momentum, which is intrinsically provided by the formulation of the Hamiltonian (1). In our derivation, in fact, we did not consider any contribution to the energy that may arise from the angular momentum carried by each mode, i.e., analogously to rotational kinetic energy of classical mechanics. Besides, we studied the input power dependence of the Hamiltonian and the total number of particles, whose conservation laws are at the basis of our theoretical model. Since we are injecting the laser beam onto the fiber axis, the conserved angular momentum of the beam remains equal to zero. Hence, it does not play a role in our study. On the other hand, whenever beams carrying orbital angular momentum are coupled to the fiber (e.g., when Gaussian beams are injected with a tilt angle with respect to the fiber axis and an off-set with respect the fiber core), a generalized equilibrium distribution can be derived. This has been recently shown in Refs. [48–50], where the conservation of the angular momentum is considered as a further constraint to be respected in the thermalization process.
In Fig. 9, we show that both the Hamiltonian and the parity remain nearly constant when varying the input peak power, for all input pulse durations. Specifically, we found that $H_{exp}\simeq 55.3$ $\mathrm{mm}^{-1}$ and $M_{exp} \simeq 0.04$, $H_{exp}\simeq 55.4$ $\mathrm{mm}^{-1}$ and $M_{exp} \simeq -0.02$, $H_{exp}\simeq 54.1$ $\mathrm{mm}^{-1}$ and $M_{exp} \simeq 0.01$, at 174 fs, 1 ps, and 7.6 ps of pulse duration, respectively. These results can be compared to the values of $H$ and $M$ that are calculated by substituting the RJ distribution which fits the experimental data in Eq. (1) and Eq. (8). We found that $H_{RJ} = 56.7$ $\mathrm{mm}^{-1}$, $H_{RJ}\simeq 57.6$ $\mathrm{mm}^{-1}$, $H_{RJ}\simeq 56.1$ $\mathrm{mm}^{-1}$ at 174 fs, 1 ps, and 7.6 ps of pulse duration, respectively. Whereas $M_{RJ} = 0$, owing to the above-mentioned symmetry of the RJ distribution.
Fig. 9. Conservation of the Hamiltonian $H$ and of parity $M$, for an input pulse duration of 174 fs (a), 1 ps (b) and 7.6 ps (c), respectively. The error bars are estimated by considering all of the reconstructions of the output beam near-field at each input power (see the MD method in the Appendix Section). The horizontal dashed lines are obtained from the RJ distribution that fits the experimental data in the BSC regime.
Finally, the number of particles is well-known to be conserved, since linear absorption is less than 1 dB/km at the wavelength of operation [51]. The number of particles conservation law would be broken if we were working in different regimes, e.g., when dealing with non-negligible linear or nonlinear absorption. Nevertheless, this would require to reach, on one side, ultraviolet or mid-infrared wavelengths (i.e. for reaching the silica bandgap or triggering phonon resonances) or, on the other, pulse peak powers that are well above the BSC threshold [52].
5. Conclusions
In this work we provided a comprehensive description, based on statistical mechanics, of the physical mechanisms behind BSC. Our approach permits a simple derivation, from fundamental principles, of the mode distribution at the thermodynamic equilibrium. We have compared the theoretical predictions of the BSC mode distributions in GRIN MMFs, with a full experimental characterization, based on an MD holographic technique. Theory and experiments were found to be in a very good quantitative agreement. Specifically, the mode distribution associated with BSC, that is experimentally observed at the output of a GRIN fiber when increasing the input power, can be well described by a RJ distribution, indicating that optical thermalization is reached. The parameters that define the RJ distribution, i.e., the temperature and the chemical potential, only depend on the spatial coupling conditions of the input beam. Whereas the equilibrium distribution parameters turn out to be independent of the input pulse duration. In the thermodynamic analogy, larger input powers mean greater numbers of particles, and a faster approach to an ergodic mixing of all states, or thermalization. The thermodynamic equilibrium is reached by ensuring the conservation of the total energy and the number of particles of the multimode system. Furthermore, we experimentally show that the symmetry with respect to the index $m$, which is related to the total angular momentum, is also a conserved quantity. Finally, our observations reveal that, as a result of thermalization, modes with radial symmetry exhibit a macroscopic occupation, at the expense of modes with more complex azimuthal characteristics. The results of our studies provide a significant contribution to the fundamental understanding of an intriguing physical process, and will be of general interest for applications of multimode nonlinear fibers in different emerging technologies, ranging from nonlinear imaging to beam delivery and high-power lasers. In conclusion, we underline that the dynamics of Hamiltonian multimode optical systems also admit a general statistical mechanics description in terms of photonic spin variables [53]. This has been experimentally demonstrated in a disordered photorefractive waveguide [54], and it can be analogously demonstrated in MMFs.
Appendix A: Mode-decomposition method
Here, we summarize the main steps of the MD method, based on digital computer holography [36]. The transverse ($x,y$) beam profile at the output of a GRIN MMF may be represented as the superposition of Laguerre-Gauss (LG) modes
where $\xi _{\ell,m}$ are the normalized complex amplitudes of each mode and $\psi _{\ell,m}$ are LG functions. The main objective of the MD is that of searching for the values of $\xi _{\ell,m}$, which are complex numbers that are defined as the scalar product between $U$ and $\psi _{\ell,m}$, i.e.,Actually, as we target the mode occupancy $n_{\ell,m}$, we do not need to know both amplitude and phase of $\xi _{\ell,m}$, but only its square modulus, as $|\xi _{\ell,m}|^{2} = n_{\ell,m}$.
In our experiments, we measure the far-field (FF) intensity $I$ of the superposition of the phase mask provided by the SLM and $U(x,y)$, appropriately magnified by the lenses $L_1$ and $L_2$ (see Fig. 3). We dub $T$ the transmission function which transforms $U$ as a consequence of $L_1-L_2$ magnification and the SLM phase encoding. Then the field profile on the FF camera is obtained by calculating the Fourier transform (given by the lens $L_3$) of the product of $T$ and $U$, so that $I$ reads as
With this substitution, it is clear that if the transmission function is chosen such that $T_{\ell,m} (x,y) = \psi _{\ell,m}^{*} (x,y)$, Eq. (12) boils down to the square modulus of Eq. (10) calculated in $k_x = k_y = 0$. Therefore, by properly choosing a set of phase patterns on the SLM, a corresponding set of images captured by the FF camera can be used for calculating all of the $n_{\ell,m}$. For a full explanation of the method, and details about the phase patterns, one can refer to [40].
A.1 Estimation of the decomposition error
In Fig. 9, we show that, when varying the input power, the total angular momentum and the Hamiltonian are a constant, within the limits of the experimental error. In our report, we fully address the latter to the numerical reconstruction of the mode composition. As a matter of fact, the experimental MD analysis consists of two parts. In the first part, an algorithm simultaneously controls the SLM and the FF camera, storing then a set of images, each of which corresponds to a given phase pattern of the SLM. The second part starts from these images for reconstructing the near-field (NF) of the beam at the output facet of the fiber, by determining the parameters $n_{\ell, m}$, as discussed above. In this second part, the choice of the center ($k_x =0, k_y = 0$) is both pivotal and non-trivial. The images, in fact, are the superposition of the output beam with a series of phase patterns. The identification of a center is not necessarily straightforward, as the $I(k_x,k_y)$ does not exhibit particular symmetries (even at the occurrence of BSC).
We verified that even a single pixel offset in the choice of the center point may lead to the loss of faithful reconstructions. Therefore, we estimate the error bars in Fig. 9 by running several times the reconstruction algorithm for a given set of acquired images. We first choose the most faithful reconstruction among the several we made, thus defining the mean value of the experimental Hamiltonian and total angular momentum. Then, we define their error as the difference from the reconstructions that we obtained by choosing a center point that has a single pixel offset with respect to the chosen one.
Appendix B: Temporal and spectral broadening effects
In the main text, we only report spatial properties of the beam at the output of GRIN MMFs. This is because BSC is a spatial phenomenon, which mainly involves spatial dynamics. As we pointed out in the main text, the BSC theoretical description can be derived for CWs. Nevertheless, when dealing with pulsed laser sources, one naturally wonders what are the temporal and spectral effects during the propagation inside a nonlinear medium. Here, we consider the case of the shortest pulses used in our experiments, i.e., that having 174 fs of duration, since temporal and spectral effects are the most pronounced in this case.
As the power values involved are high-enough for triggering self-phase modulation, a spectral broadening of the initial pulse is produced. In Fig. 10.a, we report the output spectrum for several values of $P_p$. Specifically, the latter ranges from 12 kW up to 30 kW, which is the threshold value for achieving BSC. For a clearer comparison, the $P_p$ values reported here are the same as in Fig. 5. Spectra are shown in a log scale in order to emphasize their broadening. In Fig. 10(b), we show the dependence on the input peak power of the output pulse full-wave-half-maximum bandwidth ($\Delta \lambda$), calculated from Fig. 10(a). As it can be seen, $\Delta \lambda$ reaches the value of 85 nm at the occurrence of BSC. This explains the need of the band-pass filter in the experimental setup in Fig. 3 when working with short pulses. Indeed, if the filter was missing, the near-field camera would detect an incoherent superposition of waves with different wavelengths, whose average would be seen as a Gaussian "cleaned" beam. That camera, in fact, detects light with frequencies ranging from the visible to the near-infrared. Furthermore, the SLM can encode a proper phase mask only to a monochromatic source. Thus, in the absence of the band-pass filter, MD cannot be properly carried out.
Fig. 10. a) Output spectrum at input peak powers up to 30 kW (BSC power threshold). The input pulse duration is 174 fs. b) Output bandwidth vs. input peak power.
As far as temporal aspects are concerned, group velocity dispersion, in combination with self-phase-modulation has the beneficial effect of causing about tenfold broadening of the pulse temporal duration after 3 m of fiber length. As a result, modal dispersion is not able to temporally separate the input excited modes. This is shown in Fig. 10(c), where we present numerical simulations of the temporal profile of the modes at the fiber output. In the simulation, we consider the propagation over 3 m of GRIN fiber of the first 15 modes (sorted by $q$). We use the same CC as in the experiments, with input laser pulses of 174 fs, as in Fig. 5. Only three radial modes have a non-negligible power fraction at both the fiber input and output (inset of Fig. 10(c)). These modes turn out to be temporally overlapped, which justifies the fitting of our experimental results by a CW-based theoretical model.
As a final note, it is worth to underline that the conservation of the number of particles must be verified by measuring the total output power, before the output band-pass filter. For this purpose, we used a large-band thermal power meter.
Funding
European Research Council (740355); Ministero dell’Istruzione, dell’Università e della Ricerca (R18SPB8227); Ministry of Education and Science of the Russian Federation (14.Y26.31.0017); Russian Science Foundation (21-42-00019).
Acknowledgments
We acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant No. 740355), the Italian Ministry of University and Research (R18SPB8227), and the Russian Ministry of Science and Education Grant No. 14.Y26.31.0017. S.B., E.P., D.Kh. and M.G. were supported by Russian Science Foundation (Grant No. 21-42-00019).
Disclosures
The authors declare that they have no competing financial interests.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. R.-J. Essiambre, R. W. Tkach, and R. Ryf, “Chapter 1 - fiber nonlinearity and capacity: Single-mode and multimode fibers,” in Optical Fiber Telecommunications (Sixth Edition), I. P. Kaminow, T. Li, and A. E. Willner, eds. (Academic Press, Boston, 2013), Optics and Photonics, pp. 1–43, sixth edition ed.
2. D. J. Richardson, J. M. Fini, and L. E. Nelson, “Space-division multiplexing in optical fibres,” Nat. Photonics 7(5), 354–362 (2013). [CrossRef]
3. L. G. Wright, D. N. Christodoulides, and F. W. Wise, “Spatiotemporal mode-locking in multimode fiber lasers,” Science 358(6359), 94–97 (2017). [CrossRef]
4. A. Mecozzi, C. Antonelli, and M. Shtaif, “Coupled Manakov equations in multimode fibers with strongly coupled groups of modes,” Opt. Express 20(21), 23436–23441 (2012). [CrossRef]
5. S. Mumtaz, R.-J. Essiambre, and G. P. Agrawal, “Nonlinear propagation in multimode and multicore fibers: generalization of the Manakov equations,” J. Lightwave Technol. 31(3), 398–406 (2012). [CrossRef]
6. L. G. Wright, D. N. Christodoulides, and F. W. Wise, “Controllable spatiotemporal nonlinear effects in multimode fibres,” Nat. Photonics 9(5), 306–310 (2015). [CrossRef]
7. A. Picozzi, G. Millot, and S. Wabnitz, “Nonlinear virtues of multimode fibre,” Nat. Photonics 9(5), 289–291 (2015). [CrossRef]
8. K. Krupa, A. Tonello, A. Barthélémy, T. Mansuryan, V. Couderc, G. Millot, P. Grelu, D. Modotto, S. A. Babin, and S. Wabnitz, “Multimode nonlinear fiber optics, a spatiotemporal avenue,” APL Photonics 4(11), 110901 (2019). [CrossRef]
9. F. Mangini, M. Ferraro, M. Zitelli, V. Kalashnikov, A. Niang, T. Mansuryan, F. Frezza, A. Tonello, V. Couderc, A. Aceves, and S. Wabnitz, “Rainbow Archimedean spiral emission from optical fibres,” Sci. Rep. 11(1), 13030 (2021). [CrossRef]
10. K. Krupa, A. Tonello, A. Barthélémy, V. Couderc, B. M. Shalaby, A. Bendahmane, G. Millot, and S. Wabnitz, “Observation of geometric parametric instability induced by the periodic spatial self-imaging of multimode waves,” Phys. Rev. Lett. 116(18), 183901 (2016). [CrossRef]
11. K. Krupa, A. Tonello, B. M. Shalaby, M. Fabert, A. Barthélémy, G. Millot, S. Wabnitz, and V. Couderc, “Spatial beam self-cleaning in multimode fibres,” Nat. Photonics 11(4), 237–241 (2017). [CrossRef]
12. G. Lopez-Galmiche, Z. S. Eznaveh, M. A. Eftekhar, J. A. Lopez, L. G. Wright, F. Wise, D. Christodoulides, and R. A. Correa, “Visible supercontinuum generation in a graded index multimode fiber pumped at 1064 nm,” Opt. Lett. 41(11), 2553 (2016). [CrossRef]
13. Z. Liu, L. G. Wright, D. N. Christodoulides, and F. W. Wise, “Kerr self-cleaning of femtosecond-pulsed beams in graded-index multimode fiber,” Opt. Lett. 41(16), 3675–3678 (2016). [CrossRef]
14. K. Krupa, C. Louot, V. Couderc, M. Fabert, R. Guenard, B. M. Shalaby, A. Tonello, D. Pagnoux, P. Leproux, A. Bendahmane, R. Dupiol, G. Millot, and S. Wabnitz, “Spatiotemporal characterization of supercontinuum extending from the visible to the mid-infrared in a multimode graded-index optical fiber,” Opt. Lett. 41(24), 5785–5788 (2016). [CrossRef]
15. L. G. Wright, Z. Liu, D. A. Nolan, M.-J. Li, D. N. Christodoulides, and F. W. Wise, “Self-organized instability in graded-index multimode fibres,” Nat. Photonics 10(12), 771–776 (2016). [CrossRef]
16. R. Guenard, K. Krupa, R. Dupiol, M. Fabert, A. Bendahmane, V. Kermene, A. Desfarges-Berthelemot, J. L. Auguste, A. Tonello, A. Barthélémy, G. Millot, S. Wabnitz, and V. Couderc, “Kerr self-cleaning of pulsed beam in an ytterbium doped multimode fiber,” Opt. Express 25(5), 4783–4792 (2017). [CrossRef]
17. R. Guenard, K. Krupa, R. Dupiol, M. Fabert, A. Bendahmane, V. Kermene, A. Desfarges-Berthelemot, J. L. Auguste, A. Tonello, A. Barthélémy, G. Millot, S. Wabnitz, and V. Couderc, “Nonlinear beam self-cleaning in a coupled cavity composite laser based on multimode fiber,” Opt. Express 25(19), 22219–22227 (2017). [CrossRef]
18. A. Niang, D. Modotto, A. Tonello, F. Mangini, U. Minoni, M. Zitelli, M. Fabert, M. A. Jima, O. N. Egorova, A. E. Levchenko, S. L. Semjonov, D. S. Lipatov, S. Babin, V. Couderc, and S. Wabnitz, “Spatial beam self-cleaning in tapered yb-doped GRIN multimode fiber with decelerating nonlinearity,” IEEE Photonics J. 12(2), 1–8 (2020). [CrossRef]
19. K. Krupa, A. Tonello, V. Couderc, A. Barthélémy, G. Millot, D. Modotto, and S. Wabnitz, “Spatiotemporal light-beam compression from nonlinear mode coupling,” Phys. Rev. A 97(4), 043836 (2018). [CrossRef]
20. E. Podivilov, D. Kharenko, V. Gonta, K. Krupa, O. Sidelnikov, S. Turitsyn, M. Fedoruk, S. Babin, and S. Wabnitz, “Hydrodynamic 2d turbulence and spatial beam condensation in multimode optical fibers,” Phys. Rev. Lett. 122(10), 103902 (2019). [CrossRef]
21. M. Fabert, M. Săpânan, K. Krupa, A. Tonello, Y. Leventoux, S. Février, T. Mansuryan, A. Niang, B. Wetzel, G. Millot, S. Wabnitz, and V. Couderc, “Coherent combining of self-cleaned multimode beams,” Sci. Rep. 10(1), 20481 (2020). [CrossRef]
22. Y. Leventoux, A. Parriaux, O. Sidelnikov, G. Granger, M. Jossent, L. Lavoute, D. Gaponov, M. Fabert, A. Tonello, K. Krupa, A. Desfarges-Berthelemot, V. Kermene, G. Millot, S. Février, S. Wabnitz, and V. Couderc, “Highly efficient few-mode spatial beam self-cleaning at 1.5μm,” Opt. Express 28(10), 14333–14344 (2020). [CrossRef]
23. Y. Wu, H. Pourbeyram, D. N. Christodoulides, and F. W. Wise, “Weak beam self-cleaning of femtosecond pulses in the anomalous dispersion regime,” Opt. Lett. 46(13), 3312–3315 (2021). [CrossRef]
24. M. Zitelli, M. Ferraro, F. Mangini, and S. Wabnitz, “Single-mode spatiotemporal soliton attractor in multimode GRIN fibers,” Photonics Res. 9(5), 741–748 (2021). [CrossRef]
25. Z. Mohammadzahery, M. Jandaghi, S. Alipour, S. S. Rizi, E. Hajinia, E. Aghayari, and H. Nabavi, “Nonlinear spatial reshaping of pulsed beam in a step-index few-mode optical fiber,” Opt. Express 29(7), 10716–10725 (2021). [CrossRef]
26. Y. Leventoux, G. Granger, K. Krupa, A. Tonello, G. Millot, M. Ferraro, F. Mangini, M. Zitelli, S. Wabnitz, S. Février, and V. Couderc, “3D time-domain beam mapping for studying nonlinear dynamics in multimode optical fibers,” Opt. Lett. 46(1), 66–69 (2021). [CrossRef]
27. F. O. Wu, A. U. Hassan, and D. N. Christodoulides, “Thermodynamic theory of highly multimoded nonlinear optical systems,” Nat. Photonics 13(11), 776–782 (2019). [CrossRef]
28. H. Pourbeyram, P. Sidorenko, F. O. Wu, D. N. Christodoulides, and F. W. Wise, “Optical Thermalization in Ultrashort Pulse Propagation in Multimode Fiber,” in Conf. Lasers Electro-Optics, (OSA, Washington, D.C., 2020), p. SM1P.6).
29. H. Pourbeyram, P. Sidorenko, F. Wu, L. Wright, D. Christodoulides, and F. Wise, “Direct measurement of thermalization to Rayleigh-Jeans distribution in optical beam self-cleaning,” arXiv preprint arXiv:2012.12110 (2020).
30. P. Aschieri, J. Garnier, C. Michel, V. Doya, and A. Picozzi, “Condensation and thermalization of classsical optical waves in a waveguide,” Phys. Rev. A 83(3), 033838 (2011). [CrossRef]
31. K. Baudin, A. Fusaro, K. Krupa, J. Garnier, S. Rica, G. Millot, and A. Picozzi, “Classical Rayleigh-Jeans condensation of light waves: Observation and thermodynamic characterization,” Phys. Rev. Lett. 125(24), 244101 (2020). [CrossRef]
32. E. Deliancourt, M. Fabert, A. Tonello, K. Krupa, A. Desfarges-Berthelemot, V. Kermene, G. Millot, A. Barthélémy, S. Wabnitz, and V. Couderc, “Kerr beam self-cleaning on the lp11 mode in graded-index multimode fibers,” OSA Continuum 2(4), 1089–1096 (2019). [CrossRef]
33. E. Deliancourt, M. Fabert, A. Tonello, K. Krupa, A. Desfarges-Berthelemot, V. Kermene, G. Millot, A. Barthélémy, S. Wabnitz, and V. Couderc, “Wavefront shaping for optimized many-mode kerr beam self-cleaning in graded-index multimode fiber,” Opt. Express 27(12), 17311–17321 (2019). [CrossRef]
34. L. Li, J. Leng, P. Zhou, and J. Chen, “Multimode fiber modal decomposition based on hybrid genetic global optimization algorithm,” Opt. Express 25(17), 19680–19690 (2017). [CrossRef]
35. C. Schulze, D. Naidoo, D. Flamm, O. A. Schmidt, A. Forbes, and M. Duparré, “Wavefront reconstruction by modal decomposition,” Opt. Express 20(18), 19714–19725 (2012). [CrossRef]
36. D. Flamm, D. Naidoo, C. Schulze, A. Forbes, and M. Duparré, “Mode analysis with a spatial light modulator as a correlation filter,” Opt. Lett. 37(13), 2478–2480 (2012). [CrossRef]
37. D. Schimpf, R. Barankov, and S. Ramachandran, “Cross-correlated χ2 imaging of fiber and waveguide modes,” Opt. Express 19(14), 13008–13019 (2011). [CrossRef]
38. G. Pariente, V. Gallet, A. Borot, O. Gobert, and F. Quéré, “Space–time characterization of ultra-intense femtosecond laser beams,” Nat. Photonics 10(8), 547–553 (2016). [CrossRef]
39. J. Nicholson, A. D. Yablon, S. Ramachandran, and S. Ghalmi, “Spatially and spectrally resolved imaging of modal content in large-mode-area fibers,” Opt. Express 16(10), 7233–7243 (2008). [CrossRef]
40. M. Gervaziev, I. Zhdanov, D. Kharenko, V. Gonta, V. Volosi, E. Podivilov, S. Babin, and S. Wabnitz, “Mode decomposition of multimode optical fiber beams by phase-only spatial light modulator,” Laser Phys. Lett. 18(1), 015101 (2020). [CrossRef]
41. F. Mangini, M. Ferraro, M. Zitelli, A. Niang, A. Tonello, V. Couderc, O. Sidelnikov, F. Frezza, and S. Wabnitz, “Experimental observation of self-imaging in SMF-28 optical fibers,” Opt. Express 29(8), 12625–12633 (2021). [CrossRef]
42. H. Touchette, “Equivalence and nonequivalence of ensembles: thermodynamic, macrostate, and measure levels,” J. Stat. Phys. 159(5), 987–1016 (2015). [CrossRef]
43. E. M. Lifshitz and L. P. Pitaevskii, Statistical physics: theory of the condensed state, vol. 9 (Elsevier, 2013).
44. L. Rishøj, B. Tai, P. Kristensen, and S. Ramachandran, “Soliton self-mode conversion: revisiting raman scattering of ultrashort pulses,” Optica 6(3), 304–308 (2019). [CrossRef]
45. T. Hansson, A. Tonello, T. Mansuryan, F. Mangini, M. Zitelli, M. Ferraro, A. Niang, R. Crescenzi, S. Wabnitz, and V. Couderc, “Nonlinear beam self-imaging and self-focusing dynamics in a GRIN multimode optical fiber: theory and experiments,” Opt. Express 28(16), 24005–24021 (2020). [CrossRef]
46. K. G. Makris, F. O. Wu, P. S. Jung, and D. N. Christodoulides, “Statistical mechanics of weakly nonlinear optical multimode gases,” Opt. Lett. 45(7), 1651–1654 (2020). [CrossRef]
47. K. Krupa, G. G. Castañeda, A. Tonello, A. Niang, D. S. Kharenko, M. Fabert, V. Couderc, G. Millot, U. Minoni, D. Modotto, and S. Wabnitz, “Nonlinear polarization dynamics of kerr beam self-cleaning in a graded-index multimode optical fiber,” Opt. Lett. 44(1), 171–174 (2019). [CrossRef]
48. F. O. Wu, Q. Zhong, H. Ren, P. S. Jung, M. Khajavikhan, and D. N. Christodoulides, “Thermalization of orbital angular momentum in highly multimoded nonlinear optical fibers,” in 2021 Conference on Lasers and Electro-Optics (CLEO), (IEEE, 2021), pp. 1–2.
49. E. Podivilov, F. Mangini, O. Sidelnikov, M. Ferraro, M. Gervaziev, D. Kharenko, M. Zitelli, M. Fedoruk, S. Babin, and S. Wabnitz, “Thermalization of orbital angular momentum beams in multimode optical fibers,” arXiv preprint arXiv:2112.13696 (2021).
50. F. O. Wu, Q. Zhong, H. Ren, P. S. Jung, K. G. Makris, and D. N. Christodoulides, “Thermalization of light’s orbital angular momentum in nonlinear multimode waveguide systems,” arXiv preprint arXiv:2112.10841 (2021).
51. G. P. Agrawal, “Nonlinear fiber optics,” in Nonlinear Science at the Dawn of the 21st Century, (Springer, 2000), pp. 195–211.
52. M. Ferraro, F. Mangini, M. Zitelli, A. Tonello, A. De Luca, V. Courderc, and S. Wabnitz, “Femtosecond nonlinear losses in multimode optical fibers,” Photonics Res. 9(12), 2443 (2021). [CrossRef]
53. A. Ramos, L. Fernández-Alcázar, T. Kottos, and B. Shapiro, “Optical phase transitions in photonic networks: a spin-system formulation,” Phys. Rev. X 10(3), 031024 (2020). [CrossRef]
54. D. Pierangeli, A. Tavani, F. Di Mei, A. J. Agranat, C. Conti, and E. DelRe, “Observation of replica symmetry breaking in disordered nonlinear wave propagation,” Nat. Commun. 8(1), 1501 (2017). [CrossRef]
