Spatiotemporal mode decomposition of ultrashort pulses in linear and nonlinear graded-index multimode fibers

Mario Zitelli; Vincent Couderc; Mario Ferraro; Fabio Mangini; Pedro Parra-Rivas; Yifan Sun; Stefan Wabnitz

doi:10.1364/PRJ.484271

Abstract

We develop a spatiotemporal mode decomposition technique to study the spatial and temporal mode power distribution of ultrashort pulses in long spans of graded-index multimode fiber, for different input laser conditions. We find that the beam mode power content in the dispersive pulse propagation regime can be described by the Bose–Einstein law, as a result of the process of power diffusion from linear and nonlinear mode coupling among nondegenerate mode groups. In the soliton regime, the output mode power distribution approaches the Rayleigh–Jeans law.

1. INTRODUCTION

In recent years, there has been a resurgence of interest in studying nonlinear effects in multimode optical fibers (MMFs), both for their application and fundamental science perspectives [1,2]. Among these, the surprising observation of spatial beam self-cleaning induced by the Kerr effect in graded-index (GRIN) MMFs has attracted significant attention [3,4]. Analogous, but different in nature from Raman beam clean-up [5], beam self-cleaning consists of the reshaping, as the input power is increased, of the multimode speckled transverse intensity profile at the output of the fiber into a bell-shaped beam close to the fundamental mode of the fiber, sitting on a low-intensity multimode background [3].

Spatial beam self-cleaning can be exploited to substantially improve the beam quality at the output of MMFs, a property which is of interest for high-peak power beam delivery, mode-locked fiber lasers [6,7], and high-resolution imaging systems [8,9], to cite a few. To date, several models have been proposed for describing the physical mechanism of Kerr self-cleaning, including nonlinear nonreciprocal mode coupling [3], self-organized instabilities [10], 2D hydrodynamic turbulence [11], beam condensation [12,13], and, last but not least, beam thermalization [14]. In a thermodynamic framework, the multimode beam is represented as a gas of indistinguishable particles, that populates the large, but finite number of spatial modes of the fiber. This approach is powerful, since it permits the exploitation of the conservation laws (e.g., the number of particles, or power, the linear momentum, and the orbital angular momentum), which are inherent to virtually lossless beam propagation in short lengths of fiber. By exploiting the tools of statistical mechanics, one may then describe beam self-cleaning as the irreversible evolution of the gas of photons toward its state of thermal equilibrium. This is identified by the definition of thermodynamic parameters, such as the temperature ( $T$ ) and the chemical potential ( $μ$ ) [14]. The validity of this approach has been recently confirmed by means of mode-decomposition experiments, revealing the input power dependence of the mode redistribution at the fiber output, in good qualitative agreement with theoretical predictions [15 –17].

In the thermodynamic approach, the role of nonlinearity is merely that of inducing, via four-wave mixing, the exchange of photons among nondegenerate fiber modes, in analogy with particle collisions in a gas [14]. However, it has been noted that mode energy exchanges leading to thermal equilibrium are greatly facilitated by the presence of random linear mode coupling (RMC) [18,19]. The question then arises: is nonlinearity a necessary ingredient for thermalization in highly MMFs, or can linear disorder do the job by itself? And in this case, what is the expected output mode power distribution? In Raman fiber laser experiments based on kilometer-long spans of GRIN MMF, it has been observed that the highly multimode pump beam reaches a steady-state distribution that corresponds to mode equipartition [20 –22]; whereas in a recent work, the roles of nonlinear and linear mode coupling have been systematically studied, by artificially varying the strength of linear disorder in an MMF [23]. These experiments confirmed that linear and nonlinear coupling processes may lead to different thermalization regimes, which are consistent with either the mode equipartition (when linear disorder prevails) or a Rayleigh–Jeans (RJ) distribution (when nonlinearity dominates). In Ref. [23], strong linear disorder was introduced by applying a stress to the MMF via a series of clamps. It remains to be studied what is the role of the linear disorder which is inherent to standard telecommunication GRIN fibers, owing to the presence of microbends that lead to nearest-neighbor coupling of power among adjacent nondegenerate modes [24], and its interplay with a weak fiber nonlinearity.

Sign up for Photonics Research TOC. Get the latest issue of Photonics Research delivered right to you！Sign up now

This is the purpose of this work: by applying a novel time-resolved mode decomposition technique, which permits measuring with unprecedented accuracy the energy fraction of higher-order modes, we could determine the evolution of the mode power distribution in long lengths (up to 5 km) $L$ of GRIN MMF. We explored four propagation regimes: (i) a weakly nonlinear dispersive pulse regime, where the effective interaction length $L_{e}$ was of the same order as the nonlinear length, as in the early self-cleaning experiments [3], but with $L_{e} ≪ L$ ; (ii) a pseudo-linear pulsed regime [25], where strong dispersion-induced pulse broadening rapidly imposes linear propagation; (iii) a purely linear continuous-wave (CW) regime; and (iv) a nonlinear solitonic propagation regime.

In the former three cases, linear random mode coupling prevails over weakly nonlinear mode mixing in determining the steady-state mode power distribution, which is a realistic case in most applications. Surprisingly, our observations reveal that, even in these regimes, no mode equipartition was observed in GRIN fibers. To the contrary, the long-distance mode distribution converges toward a thermalized state, that is quantitatively well described by an exponential decay versus mode energy, as described by a Bose–Einstein (BE) distribution. Note that a similar exponential mode power distribution was recently observed for the Stokes beam in a Raman fiber laser; however, in that case mode-dependent gain has a significant influence in determining the equilibrium mode distribution [22]. In the soliton case, beam condensation in the fundamental fiber mode and pulse narrowing were observed. Here the mode power distribution approaches the Rayleigh-Jeans law, albeit with a slightly larger depletion of higher-order modes. Our experiments revealed the occurrence of a sudden change of the output modal distribution, from the BE to a quasi-RJ law, when a temporal soliton is forming. This finding provides, we believe for the first time, an experimental evidence of a distribution transition, which occurs in long spans of multimode optical fiber when passing from a linear-disorder dominated and pseudo-linear dispersive propagation regime, to a soliton regime.

2. THEORY

Modes of GRIN fibers are defined in terms of radial and azimuthal numbers ( $ℓ, m$ ), and can be grouped in terms of their quantum number $q = 2 ℓ + | m | + 1$ , with $q = 1, 2, \dots, Q$ because of their equal propagation constant, $β_{q}$ . Besides the total number of power packets $N = \sum n_{q}$ , the total momentum or energy $H = \sum ϵ_{q} n_{q}$ , where $ϵ_{q} = β_{q} - β_{Q}$ , is also conserved upon beam propagation. Maximization of the entropy, constrained by the conservation of $N$ and $H$ , leads to the BE law [14,16]: ${| c_{q} |}^{2} = \frac{n_{q}}{g_{q}} = \frac{1}{\exp [- (μ + ϵ_{q}) / T] - 1},$ (1)where $g_{q} = q$ is the number of modes within group $q$ . Whenever $| (μ + ϵ_{q}) / T | ≪ 1$ , Eq. (1) can be approximated by the RJ distribution ${| c_{q} |}^{2} = - T / (μ + ϵ_{q})$ . As we shall see, for high-order modes (i.e., with small $ϵ_{q}$ ) the approximation leading to the RJ law becomes invalid, meaning that the experimental mode distribution is described by the original BE law Eq. (1). In the limit of high temperatures, the mode distribution can be approximated with the Maxwell–Boltzmann (MB, or exponential) law ${| c_{q} |}^{2} ≃ \exp [(μ + ϵ_{q}) / T]$ , i.e., the equilibrium distribution of a gas of distinguishable classical particles. In the experiments described in the next section, mode group powers $n_{q}$ were estimated from fast photodiode traces. The average individual mode power fraction was determined as ${| c_{i} |}^{2} = n_{q} / q$ , with $i = 1$ for $q = 1$ , $i = 2$ , 3 for $q = 2$ , $i = 4$ , 5, 6 for $q = 3$ , and so on. This is justified by considering that RMC randomly scrambles the power of degenerate modes; hence, a statistical equipartition of mode power within each group can be assumed. Distributions of ${| c_{i} |}^{2}$ for the individual modes will be compared with BE, RJ, and MB laws, limiting the analysis to $Q = 10$ groups comprising $N = 55$ modes.

In order to verify the accuracy of a mode decomposition, it is necessary to perform a reconstruction of the output near-field intensity pattern. When using 2D modal decomposition methods such as in Ref. [26], it is challenging to properly reconstruct the output near field of ultrashort pulses propagating over long fiber lengths, because of the temporal variation of the mode envelope phases. Such variations are due to (i) the chromatic dispersion-induced phase chirp of pulses carried by different modes; (ii) the mode dispersion-induced group delay among modes; and (iii) the laser-induced phase noise.

To overcome these difficulties, it is necessary to adopt a 3D mode field reconstruction method which accounts for all of the previously listed effects, which enables us to reconstruct both the near-field transverse intensity pattern and the temporal power traces, emerging from long spans of GRIN fibers. Our 3D field reconstruction method uses the experimental output power fractions ${| c_{i} |}^{2}$ , in order to provide a weight to the numerically propagated output spatiotemporal mode fields. We simulated the linear propagation, in the absence of fiber nonlinearity and mode coupling, of each individual mode field. This provides the estimated time-dependent output mode phase profile, which accounts for the cumulated chromatic and modal dispersion, and laser linewidth; the relative phase among different modes was randomly chosen, as expected from linear RMC. We suppose that each mode carries a Gaussian pulse with initial complex amplitude $A_{p} (0, t)$ and 3D field: $E_{i} (x, y, z = 0, t) = A_{i} (0, t) F_{i} (x, y) = \sqrt{P_{i}} \exp (i ϕ_{i}) \exp [- \frac{{(t - t_{0})}^{2}}{2 T_{0}^{2}}] F_{i} (x, y),$ (2)where $F_{i} (x, y)$ is the mode eigenfunction, $P_{i}$ is the input peak power, $T_{FWHM} = 1.665 T_{0}$ is the input pulse width, and $ϕ_{i}$ is an input random phase, modeling the effect of linear RMC during propagation. In order to simulate the CW beam propagation regime, the Gaussian pulse was replaced by a laser noise term $\exp [i \int_{- \infty}^{t} ω_{0} (t^{'}) d t^{'}]$ , where the laser linewidth noise $ω_{0} (t)$ has a zero mean and standard deviation $σ_{ω} = 2 π B$ , with $B$ the laser linewidth.

Pulse envelopes $A_{i} (z, t)$ were numerically propagated according to the different modal and chromatic dispersion, and with constant linear loss [27]. The resulting output spatiotemporal field can be written as $E_{tot} (x, y, z, t) = \sum_{i = 1}^{N} c_{i} E_{i} (x, y, z, t)$ . By integrating ${| E_{tot} |}^{2}$ versus time, one may obtain the output near field. By integrating it versus $(x, y)$ , the instantaneous power profile can also be reconstructed.

3. EXPERIMENTAL RESULTS

In our experiments, we used relatively long (up to 5 km) spans of commercial GRIN OM4 fiber, with $50 / 125 μm$ core/cladding diameter, cladding index $n_{clad} = 1.444$ at 1550 nm, and relative index difference $Δ = 0.0103$ . The chromatic dispersion was $β_{2} = - 28.8 {ps}^{2} / km$ for the fundamental mode at 1550 nm, with a dispersion slope $β_{3} = 0.153 {ps}^{3} / km$ . The corresponding self-imaging (SI) period is $Λ = π a / \sqrt{2 Δ} = 0.55 mm$ , where $a$ is the core radius. Fibers were spooled on a support with 8 cm radius of curvature; care was taken to avoid stronger bending. According to Gloge theory [28], bending losses remain negligible up to the first 10 mode groups ( $α ≃ 8 \times 10^{- 10} dB / km$ for group 10, and $α ≫ 30 dB / km$ for group 11).

Several coherent sources were alternatively used at the fiber input: a 70 fs pulse source at 1550 nm wavelength, 100 kHz repetition rate, obtained from an optical parametric amplifier fed by a femtosecond Yb laser; a 1.4 ps unchirped pulse source, obtained from the 70 fs pulses by narrow bandpass filtering; a CW distributed feedback (DFB) laser with 1 nm linewidth at 1460 nm. The input beam was attenuated, linearly polarized, and passed through a $λ / 4$ wave plate, in order to generate a circular state of polarization. The Gaussian shaped beam was then injected into GRIN fiber spans of variable lengths (from 1 m up to 5 km), with a waist $w_{0}$ of approximately 15 μm, nearly twice the fundamental mode waist. The corresponding SI-induced beam compression factor $C = 2 Λ / (π β_{0} w_{0}^{2}) = 0.262$ , where $β_{0} = 2 π n_{0} / λ$ and $n_{0}$ is the core index [29,30]; the effective beam waist for determining the nonlinear coefficient is $w_{e} = \sqrt{C} w_{0} = 7.7 μm$ [30,31], close to the fundamental mode waist. The circular state of polarization was used at the input, in order to minimize power exchanges among polarization modes. Input and output states-of-polarization were checked by measuring the power after a rotating linear polarizer, both after the input wave plate and after 5 km of GRIN fiber; the power was constant at any angle of the polarizer, with a tolerance of 5%. In some tests, a diffuser was introduced at the input, in order to increase the input high-order mode content. Alternatively, the input of the fiber was laterally shifted with respect to the beam by up to $\pm 15 μm$ , again to increase the contribution from higher-order modes.

At the fiber output, a microlens focused the near field on an InGaAs camera (Hamamatsu C12741-03); a second lens focused the beam into a real-time multiple octave spectrum analyzer (Fastlite Mozza) with a spectral detection range of 1100–5000 nm. The output pulse temporal shape was inspected 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 was also used for the characterization of the input pulses. An $M^{2}$ tester (Gentec Beamage M2) was used to characterize the output beam quality. Input and output powers were measured by a power meter with microwatt resolution.

Figure 1 shows the results of a first test, performed with unchirped 1.4 ps pulses at 1550 nm, propagating over 5 km of GRIN fiber. No diffuser was used at the input, and the beam lateral shift was limited to 4 μm. In this case, fiber chromatic dispersion broadens the pulse up to 285 ps, whereas modal dispersion separates the mode groups by up to 9 ns. The total pulse energy in the fiber was 100 pJ, which corresponds to $P_{0} = 67 W$ of input pulse peak power. By considering the power decrease due to dispersive pulse spreading and the weak uniform linear loss, the effective nonlinear interaction length is $L_{e} = 144 m$ ; the corresponding phase shift induced by the self-phase modulation is $Φ_{NL} = γ P_{0} L_{e} = 5.51 rad$ , where $γ = 2 π n_{2} / (λ A_{eff}) = 0.57 {(W km)}^{- 1}$ and $A_{eff} = π w_{e}^{2}$ . Here we used the nonlinear coefficient of silica $n_{2} = 2.7 \times 10^{- 20} m^{2} /W$ . This means that the experiment was performed with significant nonlinearity at the fiber input ( $L_{NL} ≃ L_{e}$ ); nonlinear effects can be neglected after the effective distance $L_{e} = 144 m$ , and up to the fiber length of 5 km.

Experiments with 1 m and 5 km of GRIN MMF, 1.4 ps pulses at 1550 nm, 100 pJ pulse energy, lateral shift of 4 μm. (a) Measured near field at 1 m. (b) Measured near field at 5 km. (c) Instantaneous power at 5 km (measured and reconstructed). (d) Reconstructed near field at 5 km.

Figure 1.Experiments with 1 m and 5 km of GRIN MMF, 1.4 ps pulses at 1550 nm, 100 pJ pulse energy, lateral shift of 4 μm. (a) Measured near field at 1 m. (b) Measured near field at 5 km. (c) Instantaneous power at 5 km (measured and reconstructed). (d) Reconstructed near field at 5 km.

The measured near field after 1 m of fiber, shown in Fig. 1(a), is characterized by a highly speckled intensity pattern, denoting the presence of a large number of modes. On the other hand, after 5 km of propagation [Fig. 1(b)], a much brighter (nearly bell-shaped) beam is observed. The output power measured after 5 km is illustrated in Fig. 1(c), blue curve, showing a train of temporally broadened pulses, corresponding to the propagated mode groups (up to 10 among them are visible).

A mode decomposition (MD) for the mode powers was obtained by directly measuring the total degenerate mode group power from the photodiode response, and dividing it by the number of modes in that group; in this way one readily obtains the mean power fraction of individual modes in a degenerate group ${| c_{i} |}^{2}$ , as shown in Fig. 2(a) (black circles) against the differential propagation constant of the mode $- Δ β_{i} = β_{1} - β_{i}$ ( ${mm}^{- 1}$ ).

$(a) Mean power fractions for test with 1.4 ps, 5 km, 100 pJ, shift 4 μm. (b) Test with 70 fs, 830 m, 100 pJ, shift 15 μm. (c) T and μ for the test of case (a) versus pulse energy. (d) T and μ for the test of case (b) versus lateral shift.$

Figure 2.(a) Mean power fractions for test with 1.4 ps, 5 km, 100 pJ, shift 4 μm. (b) Test with 70 fs, 830 m, 100 pJ, shift 15 μm. (c) $T$ and $μ$ for the test of case (a) versus pulse energy. (d) $T$ and $μ$ for the test of case (b) versus lateral shift.

Figure 1(d) shows the near field after 5 km of pulse propagation, reconstructed by using the 3D method, based on the measured ${| c_{i} |}^{2}$ values from Fig. 1(c), and setting random input mode phases $0 < ϕ_{i} < 2 π$ : the image cross correlation [32] between Figs. 1(b) and 1(d) is 96%. The effect of input random phases results in a smoothing of the output near field; the images cross correlation and the near-field shape are modestly affected by a change of the input random mode phase sequence. Figure 1(c) (orange curve) is the 3D reconstruction of the output optical power, showing a good agreement with the experiment.

In a second test, unchirped 70 fs pulses at 1550 nm were propagated over 830 m of GRIN fiber; no diffuser was used at the input, and the beam lateral shift was extended up to 15 μm with respect to the fiber axis. In this case the input pulse energy was limited to 100 pJ, giving a peak power of 1.3 kW. Because of the strong dispersive spreading, the effective nonlinear interaction length is only $L_{e} = 0.64 m$ , and the overall nonlinear phase shift $Φ_{NL} = γ P_{0} L_{e} = 0.49 rad$ ; hence, $L_{NL} > L_{e}$ , and the propagation regime of this test can be considered as quasi-linear, or pseudo-linear [25].

Figure 3(a) shows the measured near field after 6 m of fiber, which reflects the highly multimodal distribution that is injected at the fiber input. After 830 m of propagation, the measured near field of Fig. 3(b) shows a homogeneous, bell-shaped output beam. The experimental instantaneous power after 830 m in Fig. 3(c), blue curve, shows a 1100 ps pulse with a trailing tail. Here chromatic dispersion is responsible for pulse broadening up to 947 ps, whereas modal dispersion introduces a 0.19 ps/m delay among adjacent mode groups. Experimentally, the overall dispersion introduces a total pulse broadening of up to 1100 ps. As a result, the different pulses belonging to nondegenerate mode groups remain strongly overlapped in time at the fiber output. Yet, it is still possible to recover the mode power distribution by sampling the measured instantaneous power at the appropriate modal delays for each mode group. After dividing these values by the number of modes in each group, one obtains the mode power distribution of Fig. 2(b) (black circles). Again, the measured mode power distribution permits, in combination with the 3D simulation method and by setting random mode phases, the reconstruction of the transverse near-field intensity profile after 830 m of GRIN fiber, as shown in Fig. 3(d). The corresponding simulated temporal power profile is reported as an orange curve in Fig. 3(c). The good reconstruction in both spatial and temporal dimensions confirms the correct sampling of the mode group powers in Fig. 3(c).

Figure 3.Experiment with 6 and 830 m of GRIN MMF, 70 fs pulses at 1550 nm, 100 pJ pulse energy, and lateral input beam offset of 15 μm. (a) Measured near field at 6 m. (b) Measured near field at 830 m. (c) Instantaneous power at 830 m (measured and reconstructed). (d) Reconstructed near field at 830 m.

Figures 2(a) and 2(b) show that the modal mean power fractions can be properly fitted by exponential curves. Recall that experiments in Fig. 2(a) were conducted in a propagation regime which was initially dominated by nonlinearity, followed by a much longer propagation distance where linear disorder completely prevails; whereas experiments in Fig. 2(b) were only weakly affected by nonlinearity even initially. In Figs. 2(a) and 2(b) we compare experimental results against different theoretical distributions, specifically, the RJ law, the exponential MB law, and the BE law of Eq. (1). In both cases, the BE distribution better approximates the experimental data up to the ninth mode group; whereas the RJ law is not suitable for describing the output mode distribution, since it fails to fit the data after the fifth mode group. The MB exponential distribution also fits well the data. By considering the BE law as the appropriate distribution, we computed the modal temperature $T$ and chemical potential $μ$ ( ${mm}^{- 1}$ ) for the first experiment (1.4 ps pulse, 5 km) when pulse energy is increased from 100 to 400 pJ; see Fig. 2(c). As can be seen, $T$ and $μ$ remain nearly constant when increasing the nonlinearity, as it is expected from theory [14].

For the second experiment (70 fs over 830 m in the pseudo-linear regime), Fig. 2(b) illustrates the mean power fractions when the input beam is coupled with nominal 15 μm lateral shift with respect to the fiber axis: as can be seen, also in this case the BE law fits very well the data. Figure 2(d) shows that, when varying the input lateral shift from 0 to 15 μm, $T$ and $μ$ remain nearly unchanged up to the limit value of 15 μm; this indicates that the equilibrium output mode distribution is nearly independent of the input coupling conditions.

Finally, we underline that, at variance with the approximation of mode degeneration, in our experiments we found that the values of $Δ β_{i}$ are slightly different between the modes within the same group. This difference was accounted for when using Eq. (1) for fitting the experimental data in Fig. 2.

In order to test a case where beam propagation occurs in the strictly linear regime, we used a CW laser source at 1460 nm. Figure 4 shows examples of measured near-field intensity profiles after 1 m, 830 m, and 5 km of GRIN fiber, respectively. At the input, a beam diffuser was added; moreover, the beam center was nominally laterally offset by 15 μm with respect to the fiber axis. The CW input power was of only 1 μW, corresponding to $L_{NL} = 3 \times 10^{6} km$ . Meaning that in this case propagation occurs in the purely linear regime: RMC is the only source of disorder. As can be seen, the near field of Fig. 4(a) at 1 m is highly multimodal; for increasing distances, linear disorder leads to a progressive beam smoothing, producing the output of Fig. 4(c) after 5 km. Because of the operation in the CW regime, it was not possible to temporally resolve the mode groups in this case, in order to estimate the output mode power distribution. Still, the generation of bright beam with a near-field intensity profile similar to the case of Fig. 3(b) was observed even in this case.

Figure 4.Experiments with 1 m, 830 m, and 5 km of GRIN fiber, with input CW beam at 1460 nm, 1 μW power, an input beam diffuser, and a lateral shift of 15 μm.

Previous experiments involved either highly dispersive pulses, of the CW propagation regime. In the former cases, the influence of the temporal dimension on the final output power distribution can be simply taken into account by considering an effective nonlinear length. This is no longer the case when the input pulse powers increased up to the point of forming multimode GRIN solitons [31,33,34]. To explore the soliton regime and the resulting output mode power distribution, input 70 fs pulses with energies ranging from 0.2 to 1.5 nJ were propagated over 5 km of GRIN fiber. The input beam was coupled with 10 μm lateral shift with respect to the fiber axis. At 1.3 nJ of input pulse energy, the pulse starts forming a multimode soliton [34]. We observed that, for input pulse powers slightly below the soliton value, the output beam was spatially compressed, while its temporal duration broadened to a few hundreds of femtoseconds. Yet, a Raman-shifted walk-off soliton did not emerge; we will refer to this state as the quasi-soliton regime. As soon as the walk-off soliton energy was reached (around 1.5 nJ at 1550 nm, not shown), nearly all of the output pulse energy emerged in the fundamental fiber mode.

Figure 5 compares the dispersive pulse regime with the quasi-soliton pulse regime. In the latter, Fig. 5(e), there is no Raman-induced pulse delay, and a soliton has not yet separated from its residual dispersive waves; hence, a time-resolved modal decomposition is still possible. Figures 5(a)–5(c) show the measured near field, the output power (with sampled mode group powers), and the mean mode power fractions, respectively, in the dispersive pulse propagation regime. An output beam waist of 13.5 μm and a beam quality factor $M^{2} = 1.75$ were measured. In this case, the average mode power fractions, that we have extracted from the photodiode trace, are again properly fitted by a BE law. However, when propagating in the multimode quasi-soliton regime, Figs. 5(d)–5(f) show the emergence of a condensed beam, with a waist of 9.4 μm and a nearly diffraction-limited beam quality factor $M^{2} = 1.05$ . This shows that most of the soliton beam power is carried by the fundamental mode. Therefore, it is reasonable to sample the first modal group in correspondence with the pulse peak in the photodiode trace. Now, Fig. 5(f) reveals that the corresponding mode power distribution no longer follows a BE law, but rather an RJ distribution. Although, the MD experiments show that higher-order modes are slightly more depleted than the prediction based on the RJ law. We may thus conjecture that the transition from the dispersive to the soliton pulse regime is akin to a phase transition, which separates qualitatively different thermalization regimes.

$Experiments with 5 km of GRIN MMF, 70 fs pulses at 1550 nm, lateral shift of 10 μm. (a)–(c) Measured near field, output power, and mean modal power fraction for input pulse energy of 200 pJ. (d)–(f) Measured near field, output power, and mean modal power fraction for input pulse energy of 1.3 nJ, corresponding to the formation of a multimode soliton.$

Figure 5.Experiments with 5 km of GRIN MMF, 70 fs pulses at 1550 nm, lateral shift of 10 μm. (a)–(c) Measured near field, output power, and mean modal power fraction for input pulse energy of 200 pJ. (d)–(f) Measured near field, output power, and mean modal power fraction for input pulse energy of 1.3 nJ, corresponding to the formation of a multimode soliton.

4. DISCUSSION AND CONCLUSIONS

In this work, a time-resolved MD technique was used, in order to analyze the output mode power distributions after long spans of GRIN fibers. We considered four different propagation regimes, i.e., purely linear, pseudo-linear, dispersive nonlinear, and soliton regimes. Our 3D MD and reconstruction method accounts for phase chirping, modal delay, and random complex mode scrambling, induced by propagation on the complex amplitude of the spatiotemporal field. The accuracy of the method was confirmed by successfully reconstructing the output field, both in the spatial and in the temporal domains. Comparisons among the experimental MD based on fast photodiode traces, the 3D modal reconstructions, and fitting of the output modal content, have shown that the BE law is the most appropriate for describing long-distance multimode fiber propagation for both dispersive nonlinear and quasi-linear regimes, in the presence of random mode coupling. Similar results in terms of near-field transverse intensity profiles were obtained in a purely linear propagation regime, although no direct confirmation of the output mode power distribution was possible in this case.

It was also observed that the multimode soliton regime is characterized by a different mode power distribution, which approaches the RJ law. The BE distribution, observed in a regime where linear disorder prevails over a weak nonlinearity, reveals the presence of a high-temperature thermalization process, where the input power is diffused out of the fundamental mode into higher-order modes. To the contrary, spatiotemporal beam condensation is approached in the soliton regime, and the output mode power distribution approaches the RJ law in this case, in agreement with former studies of Kerr-driven spatial beam self-cleaning. The difference in the two observed distributions could be related to a glassy regime [35], caused by the nonlinearity affected by the linear disorder; this property will be investigated in a future work.

We underline that the exponential decay of the mode power fraction must not be confused with the classical Maxwell–Boltzmann law. Indeed, in the latter, the gas particles are supposed to be distinguishable. Whereas our experiments involve indistinguishable photons in the fiber modes, whose thermalization process occurs before the mode distinguishability provided by the modal walk-off.

Acknowledgment

Acknowledgment. We acknowledge helpful discussions with G. Steinmeyer about using either the RJ or the BE law for describing the mode power distribution of high-order fiber modes.

References

[1] K. Krupa, A. Tonello, A. Barthélémy, T. Mansuryan, V. Couderc, G. Millot, P. Grelu, D. Modotto, S. A. Babin, S. Wabnitz. Multimode nonlinear fiber optics, a spatiotemporal avenue. APL Photonics, 4, 110901(2019).

[2] I. Cristiani, C. Lacava, G. Rademacher, B. J. Puttnam, R. S. Luìs, C. Antonelli, A. Mecozzi, M. Shtaif, D. Cozzolino, D. Bacco, L. K. Oxenløwe, J. Wang, Y. Jung, D. J. Richardson, S. Ramachandran, M. Guasoni, K. Krupa, D. Kharenko, A. Tonello, S. Wabnitz, D. B. Phillips, D. Faccio, T. G. Euser, S. Xie, P. S. J. Russell, D. Dai, Y. Yu, P. Petropoulos, F. Gardes, F. Parmigiani. Roadmap on multimode photonics. J. Opt., 24, 083001(2022).

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

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

[5] K. S. Chiang. Stimulated Raman scattering in a multimode optical fiber: evolution of modes in Stokes waves. Opt. Lett., 17, 352-354(1992).

[6] L. G. Wright, D. N. Christodoulides, F. W. Wise. Spatiotemporal mode-locking in multimode fiber lasers. Science, 358, 94-97(2017).

[7] U. Tegin, B. Rahmani, E. Kakkava, D. Psaltis, C. Moser. Single-mode output by controlling the spatiotemporal nonlinearities in mode-locked femtosecond multimode fiber lasers. Adv. Photonics, 2, 056005(2020).

[8] N. O. Moussa, T. Mansuryan, C.-H. Hage, M. Fabert, K. Krupa, A. Tonello, M. Ferraro, L. Leggio, M. Zitelli, F. Mangini, A. Niang, G. Millot, M. Papi, S. Wabnitz, V. Couderc. Spatiotemporal beam self-cleaning for high-resolution nonlinear fluorescence imaging with multimode fiber. Sci. Rep., 11, 1(2021).

[9] Y. Zhang, Y. Cao, J.-X. Cheng. High-resolution photoacoustic endoscope through beam self-cleaning in a graded index fiber. Opt. Lett., 44, 3841-3844(2019).

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

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