Fig. 1. (a) Temporal solitons propagate within a singlemode fiber, forming due to the equilibrium between chromatic dispersion induced temporal broadening and nonlinearity induced pulse compression. (b) Spatial solitons manifest in Kerr nonlinear media, shaping a beam with a consistent beam waist as it travels, emerging from the interplay of spatial diffraction induced beam divergence and nonlinear self-focusing. (c) Spatiotemporal MMSs propagate in graded-index multimode fibers, where the interaction between beam diffraction and fiber confinement creates spatial multimodal beams. MMS formation is more complex, necessitating a balance not only between chromatic dispersion broadening and nonlinear pulse compression, but also a balance between nonlinear beam trapping and modal walk-off due to the disparate group velocities of the fiber modes.

Fig. 2. Overview of the relationships between the fundamental properties of MMFs and the key physical phenomena leading to the formation of an MMS.

Fig. 3. Refractive index profiles n(x,y=0) (black solid line) of SIMF in (a), (b) and GIMF in (c), (d) and their corresponding eigenmode profiles Fp(x,y=0) (p=1,6,15), modal effective indices neff,p (p=1,…,15). (e)–(h) Mode dispersion comparison between the SIMF and the GIMF with the same core radius of 50 μm at 1.45 μm, for (e) mode effective areas, (f) relative zero-order mode dispersion δβ0(p)=β0(p)−β0,SIMF(1), (g) relative group velocity Δvg,n=1/β1(p)−1/β1,GIMF(1), (h) second-order dispersion β2(p). The fiber radius for the two fibers is R=25  μm.

Fig. 4. Characteristic lengths of GIMFs and SIMFs at different wavelengths: (a) mode beating length LB(p,q); (b) modal walk-off length LW(p,q); (c) dispersion length LD(p). In (d) we show the differential mode delay DMD(p,q) for modes p and q, for an input pulse temporal duration T0=110  fs (FWHM: 184 fs). The significance of each length scale in (a)–(c) becomes apparent when the pulse propagation distance exceeds these characteristic lengths, indicating that the corresponding effect has a significant influence. DMD(p,q) detailed in (d) quantifies the temporal separation (ps) between two modes per unit of propagation distance (m).

Fig. 5. Modal effective index neff(p)(λ) in (a) and modal group index ng(p)(λ) in (b) as functions of wavelength for mode p=1,6,15 in both GIMF and SIMF.

Fig. 6. (a) Peak power of a fundamental soliton carried by a mode with index p in SIMF and GIMF, respectively, for a 105 fs pulse duration. (b) Required peak powers and energies for forming a fundamental soliton of varying duration in mode p for SIMF and GIMF when λ=1550  nm. Fiber parameters are consistent with those in Fig. 3.

Fig. 7. Comparison of the time evolution of temporal and spectral fields in the linear and soliton regimes, respectively: (a), (b) linear regime, (c), (d) MMS regime; the input pulse has a beam size rw=3 and a duration of T0=110  fs (FWHM: 184 fs), corresponding to a dispersion length LD≈0.72  m. The normalized temporal and spectral field profiles for modes 1, 6, and 15 are depicted using red, blue, and green curves, respectively, in each panel. For (a), (b) the input pulse energy is E=0.01  nJ; for (c), (d) E=0.7  nJ; here the Raman effect is neglected, and we set fR=0.

Fig. 8. MMS within a GIMF: panels (a)–(c) illustrate the normalized temporal power for modes 1, 6, and 15, respectively, at the output of a 20 m long GIMF, versus input pulse energy. The input pulse has a Gaussian shape with a duration of T0=110  fs (FWHM: 184 fs), a wavelength of 1.45 μm, and a beam width parameter rw=3. Panels (d)–(f) detail corresponding temporal positions of peak power, peak power, and pulse duration at the fiber output. Effects of higher-order dispersion (βn(p)=0,n>2), Raman (fR=0), and self-steepening (fS=0) are disregarded, in order to focus on MMS dynamics only.

Fig. 9. Output mode evolution versus input beam size rw (input modal content): the normalized temporal field powers of (a) mode 1, (b) mode 6, and (c) mode 15 at the output of a GIMF are plotted as a function of the normalized input beam size rw, when the input pulse energy is 1 nJ, and the pulse duration is T0=110  fs (FWHM: 184 fs), corresponding to a dispersion length Ld≈2  m. Panels (d)–(f) detail the corresponding temporal positions of output peak position, peak power, and pulse duration. Here we neglect the effects of higher-order dispersion (βn(p)=0,n>2), Raman scattering (fR=0), and self-steepening (fS=0).

Fig. 10. Temporal and spectral evolution of Raman MMS in a 2 m GIMF, when the input pulse energy is (a), (b) 3 nJ and (c), (d) 7 nJ. Here fR=0.18; other parameters are the same as Fig. 7.

Fig. 11. Output spectra are plotted as a function of input pulse energy for a 2 m GIMF under different input conditions: (a) rw=2, rs=0; (b) rw=2.8, rs=0.5; (c) rw=3.2, rs=1.7. The temporal and spectral field evolutions as a function of propagation distance are depicted in (d), (e), with inset panels showing the field spectrogram at specific distances. For better visualization, the dashed region in (d) is re-plotted in (f). Panels (g) and (h) illustrate the evolution of fundamental mode content and the sum of other mode contents (modes 2 to 15) with respect to the total in (f). Reprinted with permission from Ref. [39]. Copyright 2022, Optica.

Fig. 12. Bound state, or trapped, oscillatory regime for C=2/7. Adapted from Ref. [81].

Fig. 13. (a) Regions of bound and escaping solitons in the (τ,V) phase space for θ=0. (b) Comparison of simulations carried out using Eq. (39) (solid line), Eq. (65) for θ≠0 (dashed line) and for θ=0 (dotted line). Adapted from Ref. [70].

Fig. 14. The solid line corresponds to the analytical results of Kivshar [71]. The dashed line shows the estimation of Hasegawa [33]. The numerical results of Menyuk are plotted using filled circles [73]. Here C=2/3. Adapted from Ref. [72].

Fig. 15. Bifurcation diagrams for STS states as a function of E for the anomalous (in red) and normal (in blue) GVD regimes. The green line represents the CW state of the system. (a) shows the spatial width ws of the STS as a function of E; (b) shows the inverse of the temporal width w¯t=wt−1; (c) represents the STS peak intensity Ipeak. Unstable states are represented with dashed lines, while neutrally stable by solid lines. The horizontal black line corresponds to the fundamental LG0 mode. The right column shows an example of the STS solution by using three isosurfaces at different intensities (top), and a cross-section of the same state at constant t. Adapted from Ref. [20].

Fig. 16. (a) Dependence of the energy E on κ for different values of ρ in the anomalous GVD regime. Stable (unstable) branches are plotted by using solid (dashed) lines. (b) Region of existence of STSs as a function of ρ (see green shadowed area). The line limiting that area is E=Ec. Vertical dashed lines correspond to the three cases plotted in (a). Adapted from Ref. [20].

Fig. 17. Distribution of the eigenvalues associated with the dynamical system Eq. (116) in the anomalous GVD regime. Adapted from Ref. [20].

Fig. 18. Evolution with distance z of a stable STS for E=6 [see (a)–(c)]. Panel (a) shows the variation of the peak STS intensity versus the propagation distance. Panel (b) shows a close-up view of (a) for the interval z∈[950,1000]. Panel (c) shows the evolution of the STS along the interval shown in (b), obtained by plotting two isosurfaces at I1=0.5 (red) and I2=0.1 (blue). The dashed gray straight line in (a) and (b) represents the theoretical value of the STS intensity. Adapted from Ref. [20].

Fig. 19. Panels (a) and (b) show the evolution of peak intensity of stable STS with energy E for the pure quadratic and pure quartic dispersion scenarios, respectively. In (a), the green line shows the analytical values, while the blue circles and the error bars represent the average intensity values and the standard deviation for stable states. Adapted from Ref. [87].

Fig. 20. Wave collapse started from an STS for E=9. Panel (a) shows the evolution of the Ipeak, while panel (b) shows the modification of the STS with the propagation.

Fig. 21. Temporal and spectral evolutions of an optical pulse, launched into a GRIN fiber with a peak power such that N=1 using the parameter values Cf=0.5 and q=100.

Fig. 22. Temporal and spectral evolutions of a fundamental soliton, launched into a GRIN fiber with a peak power such that N=Cf using the parameter values Cf=0.5 and q=100. Corresponding intensity profiles are compared at a distance of 0 and 10LD in the bottom panel.

Fig. 23. Temporal and spectral evolutions of a fundamental soliton over a distance of 100LD using the same parameter values used for Fig. 22. Corresponding intensity profiles are compared at a distance of 0 and 100LD in the bottom panel.

Fig. 24. Temporal and spectral evolutions of a second-order soliton over a distance of 10LD using the same parameter values used for Fig. 22. Corresponding intensity profiles are compared at a distance of 0 and 10LD in the bottom panel.

Fig. 25. Temporal and spectral evolutions of a fourth-order soliton (N=4) over one dispersion length inside a GRIN fiber with Cf=1 and q=100. Higher-order effects are included using fR=0.245 and δ3=0.02. Periodic spatial focusing of the pulsed beam does not occur for Cf=1.

Fig. 26. Temporal and spectral evolution of a fourth-order soliton (N=4) over one dispersion length inside a GRIN fiber. Parameter values are identical to those used in Fig. 25 except that the value of Cf has been reduced to 0.5. Periodic spatial focusing of the pulsed beam occurs when Cf<1 because multiple modes of the GRIN fiber are excited when Cf=0.5.

Fig. 27. Space-time plot, showing intensity of a light bullet with normalized energy E=6.

Fig. 28. Simulated evolution of peak amplitude, spatial beam width, and temporal pulse duration, respectively, for a light bullet with energy E=6.

Fig. 29. Evolution of an input hyperbolic-secant waveform A→(0,t)=A0 sech(t/t0)e^1, which in the absence of mode dispersion would result in a fundamental soliton, in a fiber supporting 2N=4 strongly coupled modes for the displayed values of the mode-dispersion coefficient κ¯MD. Time and propagation distance are normalized to the input-waveform width t0 and the fiber dispersion length LD=t02/|β¯2|, respectively.

Fig. 30. Experimental study of MMS versus input pulse energy. (a) Output spectra. (b) Measured output pulse temporal duration. (c) Measured peak wavelength. Reprinted with permission from Ref. [134]. Copyright 2013, Nature Group.

Fig. 31. Experimental demonstration that MMSs require more energy than single-mode solitons, for a given temporal duration. In the vertical axis, the figure reports experimental and simulation results for the slope of the linear relation between pulse energy and inverse temporal duration, as a function of the average size of the beam waist Rg (the waist of the fundamental mode is R0=6.5  μm). All data points are higher than the curve for single-mode solitons [solid black curve shows M1=ET0 as given by Eq. (23)]; the green point indicates the case of quasi-single-mode soliton observed in Ref. [134]. Reprinted with permission from Ref. [62]. Copyright 2015, Optical Society of America.

Fig. 32. Soliton pulse duration versus wavelength, compared with measured (green diamonds), simulated (empty red circles and blue squares, respectively), and with Eq. (208) (solid line). Reprinted with permission from Ref. [37]. Copyright 2021, Nature Group.

Fig. 33. Experimental study of high-order soliton fission in GRIN fiber. Left, autocorrelation of the output pulse versus input energy; middle, corresponding output beam profiles; right, output spectrum. Reprinted with permission from Ref. [62]. Copyright 2015, Optical Society of America.

Fig. 34. Experimental data (empty squares and empty circles for first and second Raman solitons, respectively) versus numerical simulations (solid curves) for the output soliton time width versus its energy; black crosses and black dashed line are obtained from Eq. (24) by using experimental or simulation soliton parameters, respectively. Reprinted with permission from Ref. [36]. Copyright 2020, Optical Society of America.

Fig. 35. Experimental dependence of MMS beam area versus soliton energy. Reprinted with permission from Ref. [136]. Copyright 2016, Optical Society of America.

Fig. 36. Dependence of MMS temporal width on input pulse energy for GIMF (red curve), SIMF with axial (black curve) or non-axial (blue curve) input beam. Insets show examples of autocorrelation traces. Reprinted with permission from Ref. [43]. Copyright 2019, Optical Society of America.

Fig. 37. (a) Total (black curve) and bandpass filtered (red curve) output spectrum from a step-index MMF; output spatial intensity profile before (b) and after (c) filtering; (d) spatially integrated pulse measurement; (e) simulated half width at half maximum of the spatial correlation function for randomly generated beam patterns from the last six modes (red) and the last ten modes (green), with highest-order mode indicated on x-axis; purple dashed line: half-width of the experimental spatial correlation function. Reprinted with permission from Ref. [42]. Copyright 2023, American Institute of Physics.

Fig. 38. Cutback experiments. (a) Spectra as a function of fiber length with input pulse at 1045 nm in the LP0,19 mode. (b)–(f) Output mode images from different bandpass filters. Reprinted with permission from Ref. [139]. Copyright 2019, Optical Society of America.

Fig. 39. Experimental study of MMS attractor in GIMF. Here we show the output beam waist after 1 km of GRIN fiber, versus input peak power, and 0°, 2.3°, and 4.6° input beam tilt angles; the dashed horizontal line indicates the fundamental beam waist. Insets illustrate the output beam shapes at selected powers. Reprinted with permission from Ref. [38]. Copyright 2021, Optical Society of America.

Fig. 40. Soliton collision experiments. (a) Evolution of output spectrum versus input pulse energy for a 10 m GRIN fiber. (b) Output spectrum, (c) autocorrelation, and (b) soliton temporal separation versus input energy for a 2 m GRIN fiber. Reprinted with permission from Ref. [39]. Copyright 2022, Optical Society of America.

Fig. 41. Simulated (using the MM GNLSE or the 1D NLSE model, respectively) and experimental generation of dispersive wave peaks in multimode GRIN fiber. Reprinted with permission from Ref. [41]. Copyright 2015, American Physical Society.

Fig. 42. Temporal evolution of soliton peak power in different fiber structures (in the insets, glass is depicted in gray, while white indicates empty spaces). Reprinted with permission from Ref. [147]. Copyright 2024, Elsevier.

Fig. 43. Decomposition of a Gaussian beam into the first N=30 modes of a GIMF, with panels (a)–(d) showing mode coefficients |Cp|2 against the beam lateral offset rs for beam widths rw=1,2,3,4. Radial modes (p=1,6,15,28) are predominantly excited at precise center alignment (rs=0). The blue lines in (a)–(d) indicate the total power fraction ∑p=1N|Cp|2 as a function of rs. Values below one result from the finite number of modes considered.