In the realm of nonlinear optics, the phenomenon of temporal solitons has been a subject of fascination for over five decades [1,2], ever since the first prediction of temporal solitons in 1973 [3] and their experimental observation in 1980 [4]. These self-sustained wave packets exhibit remarkable properties, propagating over long distances without distortion due to a balance between chromatic dispersion and nonlinearity [see Fig. 1(a)]. Extensive research efforts have been dedicated to unraveling the mysteries of solitons [5] and harnessing them for various applications, spanning from telecommunications [6] to supercontinuum sources [7], ultrafast mode-locked fiber laser sources [8], and cavity soliton frequency combs [9]. Similarly, spatial optical solitons, akin to their temporal counterparts, are self-trapped light beams capable of maintaining their transverse shape while propagating through nonlinear media, avoiding the divergence associated with freely diffracting beams [see Fig. 1(b)]. The presence of light alters the properties of nonlinear materials, such as refractive index, resembling the intensity profile of the beam and forming an optical lens effect. The concept of self-trapped light beams was initially introduced by Askaryan in 1962 [10], followed by comprehensive investigations in both 1D and 2D diffractive systems [11–14].

While traditional temporal or spatial solitons have primarily been explored in one or two dimensions, recent advancements have unveiled a fascinating realm of solitons involving both spatial and temporal degrees of freedom, known as spatiotemporal solitons [5,15,16], also called light bullets [17]. These solitons represent multidimensional states, as they necessitate the balance among diffraction, dispersion, and nonlinearity across spatial and temporal dimensions in the course of their propagation. This complexity presents new challenges as well as interesting opportunities in nonlinear dynamics. One significant challenge is that these multidimensional states are highly susceptible to perturbations [16]. A notable example of this fragility is the spatiotemporal wave collapse, which is commonly observed in multidimensional solitons within optical materials exhibiting Kerr nonlinearity [16,18–20]. Nevertheless, spatiotemporal solitons have been studied in various configurations [20–31]. Among them, multimode fibers (MMFs) provide a convenient test-bed [32] and have garnered considerable interest in the study of such solitons, owing to their cost-effectiveness and ease of alignment. The interaction of diffraction and spatial confinement in a multimode fiber leads to the formation of multiple co-propagating spatial modes, each with distinct dispersive properties, such as different phase and group velocities, second- and higher-order dispersions. The dispersion of these spatial modes can be balanced by Kerr nonlinearity, causing the modes to become locked together in time. As a result, the multimode field propagates as a single, shape-preserving “ball” of energy in both space and time. These equilibrium states are known as multimode solitons (MMSs) [see Fig. 1(c)].

There are two major motivations for studying MMSs: 1) the intrinsic scientific interest in propagation-invariant waves, and 2) the potential to interpret complex phenomena, especially in more than one dimension, through soliton dynamics. Interpreting instabilities, pulse formation, supercontinuum generation, and other phenomena in one dimension using the soliton concept in the frame of the 1D NLSE has proven to be highly beneficial. Similarly, MMSs could provide a foundation for understanding multidimensional nonlinear phenomena. While the community has demonstrated several instances where this approach is effective, further research is needed before we can definitively state that “a field can be described as a superposition of solitons and dispersive waves,” as is often done for one-dimensional systems.

Following early work [27,28,33–35], MMSs have gained increasing attention, particularly in graded-index (GRIN) multimode fibers (GIMFs), which minimize intermodal dispersion and reduce modal walk-off time, thus facilitating the formation of an MMS. Many experiments involving multimode solitons in MMFs have explored phenomena such as soliton fission [36], walk-off solitons [37], singlemode spatiotemporal soliton attractors [38], collisions [39], dispersive radiation sideband generation [40,41], and MMS formation in multimode step-index fibers (SIMFs) [42,43]. The development of MMS in conservative systems has spurred interest in the study of the emergence of MMS in dissipative systems, such as spatiotemporal mode-locked lasers [44–46] and different multimode cavity configurations [47,48]. Furthermore, GIMFs have emerged as valuable testbeds for highly multimode systems [21], enabling the observation of a variety of phenomena such as beam self-cleaning [49], thermalization of optical waves [50–52], and supercontinuum generation [53].

This review paper focuses on recent theoretical and experimental advancements in MMSs in conservative systems, with a primary focus on MMFs. We begin by introducing the foundational concepts of MMF and MMS, elucidating the phenomena associated with MMS in relation to the MMF properties in Section 2. Many of these phenomena, such as modal walk-off and MMS formation, can be effectively modeled by the generalized coupled mode equation, with simple examples provided in Section 3. Next, we review analytic approaches to MMS using the two-coupled mode equations in Section 4, and discuss extensions to the 3D+1 spatiotemporal nonlinear Schrödinger equation (NLSE) in Section 5. To enhance understanding, the comprehensive 3D+1 model is simplified to a reduced 1D+1 NLSE representation in Section 6, which is highly efficient for studying the self-imaging-effect-induced phenomena in GIMFs. In Section 7, we introduce an innovative Gaussian quadrature approach tailored for numerical simulations of the generalized NLSE, potentially improving computational efficiency. In Section 8, we examine the propagation regime in space-division multiplexed (SDM) optical communication systems designed with MMFs, discussing the model using a generalized Jones formalism. This model highlights the impact of linear coupling effects, primarily due to manufacturing imperfections and environmental perturbations. Finally, in Section 9, we focus on the relevant experimental studies involving the MMS phenomena.

In essence, this review endeavors to offer a comprehensive overview of MMS, from MMF properties, their theoretical analysis, innovative numerical approaches, the impact of random linear coupling, and practical experiments. By exploring these domains, the paper aims to enrich understanding and acknowledgment of the significance of MMS, thus driving progress in MMS research and technology.

In this section, we first focus on the spatial eigenmodes of two types of MMFs: the step-index multimode fiber (SIMF) and the GIMF. We discuss their key properties related to the MMS, setting the stage for understanding their impact on soliton formation. Subsequently, we introduce the generalized multimode NLSE (GMMNLSE), which is relevant for describing a broad variety of MMFs with identifiable eigenmodes. We provide concise examples, ranging from the propagation of linear multimode pulse fields to the formation of MMS, the Raman-induced MMS self-frequency shift, MMS fissions, and MMS collisions. These examples serve to effectively illustrate the discussed concepts, thereby laying a solid foundation for analytical and experimental studies.

MMFs are capable of supporting many guided modes simultaneously. To gain a comprehensive understanding of MMFs, several key concepts must be grasped, as depicted in Fig. 2. •Eigenmodes. Eigenmodes represent distinct transverse profiles of the spatial modes that can propagate within an MMF, illustrating how optical fields distribute across the fiber’s cross-section. Calculation of the modal areas and their overlap is necessary to define the modal nonlinear coefficients and the strength of nonlinear mode coupling.•Eigenvalues (neff). The modal eigenvalues, denoted as neff and representing a mode’s effective refractive index, quantify the increase in the propagation constant β=neffk0 at a specific wavelength, when compared with the vacuum wave number k0. Since neff varies with frequency, for each mode of index p the different derivatives βn(p)(ω)=∂nβ(p)∂ωn with respect to frequency (or wavelength) directly affect the properties of guided modes such as modal pulse broadening, modal group velocity, and mode beating.•χ(3)nonlinearity. The nonlinearity of an MMF results from the complex third-order optical susceptibility of the fiber’s material (typically silica). There are two main categories of cubic nonlinear effects in fibers according to their response times: the virtually instantaneous Kerr effect, and the delayed nonlinear Raman and Brillouin scattering. As the light intensity increases, the Kerr effect may compensate for both intramodal pulse broadening and intermodal walk-off, leading to the MMS formation.

Subsequent subsections will discuss each of these concepts in more detail, providing a comprehensive understanding of MMF modes and their significance in MMS formation.

Light travels longitudinally in fibers, while remaining confined transversely to the core through total internal reflection resulting from the refractive index difference between the core and the cladding. This confinement establishes boundary conditions, leading to specific guided propagation modes that depend on the fiber’s transverse refractive index profile. The diversity of index profiles gives rise to various types of MMF, each characterized by spatial modes with distinct transverse beam profiles and eigenvalues. This subsection introduces the eigenmodes by using the specific examples of the SIMF and GIMF with identical core sizes.

In the case of a SIMF, the fiber’s core has a larger refractive index nco and is surrounded by a cladding with the index ncl such that Δn=nco−ncl>0. The SIMF’s transverse refractive index profile is given by n(x,y,λ)={ncl(λ)+Δn,ifx2+y2<R2,ncl(λ),otherwise,(1)where R denotes the core’s radius. The refractive index of the fiber not only varies with the wavelength λ but also with the spatial coordinates (x,y). In the context of silica fibers, the wavelength dependence of ncl(λ) is commonly characterized using a Sellmeier equation [54], which can provide the refractive index changes over a wide wavelength range. We use the following values for the SIMF parameters: the core’s radius R=25μm and a refractive index difference Δn=0.0137 compared to the cladding. The index profile n(x,y=0,λ=1.45μm) is shown by a black curve in Fig. 3(a).

Based on the fiber’s index profile, we calculate the first 15 eigenmodes F(x,y,λs)m,l at λs=1.45μm by employing the semivectorial finite-difference method [55]. These modes F(x,y,λs)m,l, designated as linearly polarized (LP) with indices m (azimuthal) and l (radial), are displayed on the top in Fig. 3. Degenerate modes sharing the same m,l values are denoted by a, b. To simplify references to these modes, we assign a singlemode index p. For instance, F11a and F11b correspond to modes p=2 and p=3, respectively.

To illustrate the relationship among eigenmodes, eigenvalues, and refractive index distribution, we plot the mode profiles Fp(x,y=0) in Fig. 3(a) for three distinct radial modes (p=1,6,15), each shifted vertically by its corresponding mode eigenvalue, using blue and green curves at wavelengths of 1.55 μm and 1.70 μm, respectively. The corresponding modal eigenvalues at three wavelengths are shown in Fig. 3(b). These modes primarily reside within the core of the fiber. Remarkably, higher-order modes, or those with longer wavelengths, have a smaller effective index, indicating reduced confinement for them. As a result, the number of supported modes varies with wavelength. While MMFs can accommodate more modes at shorter wavelengths, they may support only a single mode at much longer wavelengths, thereby classifying them as singlemode fibers when the wavelength surpasses a certain threshold value (cutoff wavelength).

The shape of the refractive index profile of an MMF influences its guided modes. This is illustrated in Fig. 3(c) using the example of a GIMF with a core radius R=25μm. The GIMF’s refractive index follows a parabolic distribution: n(x,y,λ)={ncl(λ)+Δn(1−x2+y2R2),ifx2+y2<R2ncl(λ),otherwise.(2)

This distribution leads to eigenmodes resembling Laguerre-Gaussian (LG) modes. The refractive index n(x,y=0,λ=1.45μm) and LG mode profiles Fp(x,y=0) are depicted in Fig. 3(c) for modes p=1,6,15 (the first three radial modes: LG01, LG02, LG03) at wavelengths of 1.45 μm and 1.70 μm using red and green curves. While mode profiles for GIMF and SIMF fibers exhibit significant similarities, it is worth noting that the mode areas of a GIMF are considerably smaller than those of a SIMF. To rigorously characterize the influence of mode profiles, it is essential to define the effective areas of the MMF modes.

The mode’s area significantly influences MMF dynamics in the nonlinear (high-power) regimes. This is due to the fact that, for an identical optical power, modes with smaller spatial dimensions lead to elevated field intensities (optical power per unit area), thus inducing more pronounced refractive index changes through the Kerr effect. In the presence of nonlinear effects, the mode’s effective area (MEA) is introduced as follows [56]: Aeff(p)=(∬|Fp(x,y)|2dxdy)2∬|Fp(x,y)|4dxdy,(3)where the mode profiles are both normalized and mutually orthogonal: ∬−∞∞Fn(x,y)Fm(x,y)dxdy=δnm.(4)

The modal effective area (MEA) is not enough to quantify the strength of nonlinear effects in an MMF because the nonlinear interaction among the modes also depends on the wavelength of incident light and fiber material. Therefore, one needs to introduce the modal nonlinear coefficient (MNC), which can be calculated as γ(p)=2πn2λAeff(p),(5)where n2 is the nonlinear refractive index of the MMF that produces a change in the refractive index, Δn=n2I, which depends linearly on the local intensity I of light and has the dimensions of m2/W. Different effective mode areas in an MMF lead to different MNCs for its various modes.

The modal nonlinear length (MNL) is the characteristic length scale over which significant nonlinear optical effects occur in an optical fiber. It provides a measure of the distance along the fiber at which a light pulse is significantly affected by nonlinear interactions. The nonlinear length of a mode is defined as [56] LN(p)=[γ(p)P0(p)]−1,(6)where P0(p) is the input peak power of mode p. This form indicates that a higher input power or a larger MNC leads to a shorter nonlinear length.

It is worth noting that the nonlinear length of a mode is proportional to its MEA, indicating that the MEA is a key parameter for each mode and serves as a critical parameter for distinguishing nonlinear performance across various modes of an MMF.

Optical pulses propagating inside an MMF with high peak powers may be affected by several nonlinear effects. We briefly introduce the effects associated with the third-order nonlinear optical susceptibility (χ(3)) for which the nonlinear polarization is proportional to the third power of the electric field. Effects associated with the χ(3) nonlinearity can be grouped into two main categories according to their response time, which can be either very short or delayed in time (typically <1ps). The former encompasses the well-known Kerr effect, which produces a nearly instantaneous change in the refractive index proportional to pulse intensity. This change alters the propagation constants of all excited modes and gives rise to nonlinear phenomena such as self-phase modulation (SPM), cross-phase modulation (XPM), and four-wave mixing (FWM).

A delayed nonlinear response implies that the nonlinear polarization not only depends on the electric field at the current time, but also on the field values at earlier times. This delay is attributed to light-induced molecular vibrations involving either optical phonons (Raman scattering) or acoustic phonons (Brillouin scattering). We refer to Ref. [56] for further details. Only Raman scattering is relevant for short pulses that are of interest in this review. As a result of Raman scattering, a longer-wavelength pulse, known as the Stokes pulse, undergoes amplification in the presence of a shorter-wavelength (pump) pulse. In silica fibers, maximum gain occurs for a relative frequency shift of about 13 THz (or 100 nm) from the pump. As the intensity of the Stokes pulse increases, it can serve as the pump and generate another pulse at a longer wavelength through a cascaded Raman process, thus potentially producing multiple pulses from a single pump pulse.

The Raman gain affects ultrashort pulses (width <1ps) in a different way through a process called intrapulse Raman scattering [56]. The spectrum of such pulses is so wide that its low-frequency components can be amplified by the high-frequency components of the same pulse. Consequently, the entire spectrum of a short pulse continuously shifts towards longer wavelengths as it propagates along the fiber. This effect is particularly strong for solitons and is often referred to as the soliton self-frequency shift (SSFS). As we discuss later, it also occurs for an MMS forming inside an MMF.

Let us introduce modal dispersion, a key property of MMFs, which is associated with the eigenvalues corresponding to different guided modes. These eigenvalues are crucial, as they directly influence the propagation of optical pulses through different guided modes. We initiate our discussion by focusing on the frequency dependence of eigenvalues of the modes, and obtain the associated modal dispersion properties. By considering different orders of modal dispersion for both GIMFs and SIMFs, we are able to find several specific properties of GIMFs. GIMFs exhibit a mode-beating distance that is nearly wavelength independent, a reduced chromatic-dispersion-induced pulse broadening, and a minimal intermodal walk-off. Such attributes render GIMFs a fertile ground for observing many intriguing phenomena, including periodic self-imaging [57], formation of multimode solitons [37,38], and collisions of such solitons [39].

The eigenvalues of MMF modes are the effective modal indices np(λ). They quantify the change in a mode’s propagation constant induced by the fiber’s refractive index profile. As shown in Figs. 3(b) and 3(d), np(λ) for SIMFs and GIMFs decreases with the mode order p and wavelength λ, indicating that a weaker waveguide confinement occurs for HOMs or at longer wavelengths.

To discuss modal dispersion, we consider the frequency dependence of the propagation constant β(p)(ω) for mode p: β(p)(ω)=ωcnp(ω).(7)

Different derivative orders of β(p) with respect to frequency ω provide us with several dispersion parameters, βn(p)(ω)=∂nβ(p)∂ωn,(8)which play a crucial role when short pulses propagate inside an MMF. Utilizing these coefficients, we can elucidate the influence of dispersion on various physical phenomena in MMFs.

The propagation constant (zeroth-order dispersion) β0(p)(ω) in Eq. (8) quantifies the rate of phase change for mode p with propagation. Figure 3(f) compares for SIMF and GIMF the phase difference, δβ0(p)=β0(p)−β0,SIMF(1), between the propagation constant of mode p and that of the fundamental mode of the SIMF. In the case of a GIMF, β0(p) exhibits a distinct step-wise decrease as the mode order p increases across different degenerate mode groups.

Changes in the modal propagation constant β0(p) with p play a crucial role in MMFs and lead to phenomena such as multimode interference or mode beating. This beating rate stems from the difference of β0 for any two modes and can be used to define the mode-beating length (MBL) as LB(p,q)=2π/|β0(p)−β0(q)|.(9)

Mode beating is particularly pronounced in GIMFs, where the phase-rate differences between adjacent mode groups are nearly identical. For example, as shown in Fig. 3(f), the modes 1, 6, and 15 satisfy the relation β0(1)−β0(6)=β0(6)−β0(15). This feature enhances the mode-beating process, which involves multiple high-order modes. Whenever an optical beam excites multiple modes of a GIMF, the total intensity, obtained by adding the contributions of all modes [see Eq. (25)], exhibits a periodic evolution (the self-imaging effect). The period of mode beating is determined by the constant difference in propagation constants between adjacent groups of modes.

Interestingly, for femtosecond pulses, the intensity at points of maximum beam compression may become so high that it induces multiphoton absorption or ionization [58]. In turn, a visible photoluminescence is generated, which scatters out of the fiber as a series of bright spots, providing an experimental verification of the self-imaging phenomenon [57]. The self-imaging period does not depend on power and is equal to the beating length of different radial modes: ZSI=2π/|β0(6)−β0(1)|.(10)

By using the data plotted in Fig. 3(f), we calculate the self-imaging distance to be ZSI≈2π/(11×10−3m−1)≈0.58mm, which closely matches the experimentally observed values [57].

A remarkable property of a GIMF is its near-constant MBL for different wavelengths. To elucidate this feature, Fig. 4(a) shows the wavelength dependence of LB(1,6) for modes p=1 and p=6 for both the GIMF and the SIMF. Remarkably, as the wavelength shifts from 1.1 μm to 1.8 μm, LB in SIMF varies from 2.8 mm to 1.7 mm but its value for the GIMF remains fixed at 0.58 mm.

From a different perspective, the self-imaging period ZSI in a GIMF depends on the core radius R and the refractive index difference as ZSI≈πR/Δ [57], where Δ=(nco2−ncl2)/nco2 is the relative index difference between the core nco and the cladding ncl. Typically, the refractive index varies with wavelength, which is a common characteristic of optical materials. However, if Δ undergoes minimal changes across a range of wavelengths, the self-imaging period remains nearly constant. This feature implies that the self-imaging phenomenon in GIMFs exhibits a remarkable wavelength stability.

When non-radial modes are excited in a GIMF by launching a beam off-set from the core’s center, the beam’s shape exhibits transverse oscillations with distance. This process doubles the spatial period of self-imaging, owing to the beating between radial and non-radial modes. The resulting self-imaging period results from the propagation constant difference between two such modes and is given by ZSI′=2π/(β0(2)−β0(1))=1.16mm [see the blue curve in Fig. 4(a)].

The evolution of periodic mode beating leads to collective intensity oscillations, which in turn facilitates the emergence of spectral sidebands. This phenomenon is associated with resonant dispersive radiation in regions of anomalous dispersion [41,59] and with a so-called geometric parametric instability (GPI) in the region of normal dispersion [60]. The frequencies of the GPI-induced resonant spectral peaks adjacent to the pump frequency follow the relation fh=±h2π2πZSIβ2(1),h=1,2,3,….(11)

The group velocity of a multimode optical pulse depends not only on the optical frequency (chromatic dispersion) but also on its constituent modes. For mode p at wavelength λ, the first-order term β1(p)(λ) in Eq. (8) is inversely proportional to the group velocity vg(p)(λ)=1/β1(p)(λ).(12)

Figure 3(g) shows the relative group velocity Δvg(p)=vg(p)−vg,GIMF(p) between a mode of order p of either a GIMF or a SIMF, and the fundamental mode of a GIMF at λ=1450nm. We may observe that, when compared to SIMF, modes of GIMF propagate faster and exhibit smaller group velocity differences. Additionally, for both fibers the group velocity of LOMs is larger than that of HOMs.

It may appear surprising that LOMs have higher group velocities than HOMs, since the effective index decreases as the mode number grows larger [see Figs. 3(b) and 3(d)]. Therefore, the phase velocity of LOMs is smaller than that of HOMs. However, the group velocity has an inverse dependence on the modal group index, which includes the effects of dispersion and is defined as ng(p)(ω)=neff(p)+ωddωneff(p).(13)

For a more detailed comparison, we illustrate in Fig. 5 the wavelength dependence of both the modal effective index neff(p) and the modal group index ng(p) for modes p=1,6,15 for both GIMF and SIMF. As can be seen, contrary to the effective index, the modal group index grows larger with mode number for both types of fibers. Moreover, the modal effective index decreases monotonically with increasing wavelength. In contrast, the modal group index reaches its minimum value near the zero-dispersion wavelength, where the largest group velocity is reached. In regions of anomalous dispersion, an increasing neff for mode p means a decrease in the group velocity of that mode as its wavelength grows larger.

The differential mode delay (DMD) quantifies the temporal delay between modes per unit distance; it is defined as DMD(p,q)=β1(p)−β1(q).(14)

In Fig. 4(d) we show the wavelength dependence of DMD(p,q)(λ), for the mode pairs (1, 6) and (1, 15) for both GIMF and SIMF. For both types of fiber, the DMD increases with wavelength. In addition, an increase in the relative mode number difference leads to larger values of DMD(p,q)(λ). Consequently, the formation of MMS is much more difficult for modes with large DMDs, when compared with MMS composed of adjacent modes.

For GIMFs, the DMD values of <1ps are an order of magnitude smaller than those for SIMFs. For example, at a wavelength of 1.45 μm the temporal delay associated with modal dispersion is DMD(1,6)≈0.22ps/m for a GIMF but approximately 2.4 ps/m for a SIMF. Utilizing DMD(p,q)(λ) allows for a rough estimate of the propagation delay between modes in the linear regime. For instance, when a pulse travels 5 km inside a GIMF, its relative delay is approximately 1.1 ns for two neighboring modes. This estimate agrees with experimental findings, as found in Ref. [61].

The walk-off distance for pulses carried by different modes depends not only on the relative group velocities of these modes but also on pulse duration. Longer pulses have less pronounced walk-off compared to shorter pulses. For a pulse of duration T0, the walk-off distance between modes p and q is defined as LW(p,q)=T0/|β1(p)−β1(q)|.(15)

Here LW(p,q) stands for the distance over which two pulses separate in time by an amount equal to their temporal width. A longer walk-off length implies more interaction time due to their temporal overlap during propagation, while a shorter length suggests significant temporal separation over a shorter distance.

Low-order modes (LOMs) travel faster compared to high-order modes (HOMs), e.g., β1(1)(λ1)<β1(6)(λ1) at any given wavelength. Whenever an MMS is formed by trapping different modes via the Kerr nonlinearity, the resulting group velocity of the MMS will be intrinsically linked to its modal composition. Therefore, one may expect that an MMS composed of HOMs will propagate more slowly than an MMS composed of LOMs.

The binding mechanism for different modes composing an MMS is that the Kerr nonlinearity induces a spectral blue shift for HOMs and a red shift for LOMs. This is due to the fact that modal group velocities decrease as wavelength increases in the case of anomalous dispersion, i.e., β1(p)(λ1)<β1(p)(λ2) for λ1>λ2. Consequently, the red-shifted LOMs decelerate while the blue-shifted HOMs accelerate, gradually approaching an equilibrium point where group velocities become balanced across all modes. An illustrative example of this phenomenon can be found in Fig. 7, as detailed in Section 2.B.1.

As previously mentioned, the Raman effect induces a progressive red shift of the soliton, which decreases its group velocity. In GIMFs, the Raman effect also contributes to the MMS formation, since it reduces the speed mismatch between the fast LOMs (which have higher nonlinear coefficients and undergo the largest SSFS) and the slow HOMs. On the other hand, when two multimode Raman solitons with distinct modal compositions are generated during the fission of a single MMS, their different group velocities may lead to collisions between them as outlined in Ref. [39]. This scenario is elucidated by the example presented in Fig. 11(d) in Section 2.B.4.

Chromatic or group-velocity dispersion (GVD) characterizes the dispersive broadening of optical pulses within each mode of an MMF. As a consequence of material and waveguide dispersion, different frequency components within the pulse in a given mode travel at different speeds, resulting in pulse spreading in time upon its propagation. The GVD is governed by the second-order term β2(p)(ω) in Eq. (8). In standard silica optical fibers, there is a zero-dispersion wavelength (ZDW), close to 1.3 μm, such that β2(p)(ω)=0. Consequently, optical pulses propagating near this wavelength are little affected by GVD, but the impact of higher-order dispersions on the pulses may become significant.

For a transform-limited (i.e., unchirped) Gaussian pulse with initial temporal width T0, propagating in either the normal or anomalous dispersion regime, its temporal width increases with distance z as [56] T(z)=T0[1+(z|β2(p)|/T02)2]1/2. It is useful to introduce a characteristic dispersion length for mode p as LD(p)=T02/|β2(p)|.(16)

Consequently, for a pulse carried by mode p the pulse broadening ratio, T/T0=[1+z/LD(p)]1/2, becomes significant as it travels a distance greater than LD.

In Fig. 4(c), we plot LD(p)(λ) as a function of wavelength for both GIMF and SIMF. As can be seen, for wavelengths around the ZDW of 1.26 μm, LD(p)(λ) tends towards infinity because β2(p)(ω)≈0. The ZDW corresponds to an inflection point for β1(p)(ω), and optical pulses propagate at their fastest speed at this wavelength when compared with other wavelengths. In the normal dispersion regime (λ<1.3μm), the group velocity increases with wavelength; in the anomalous dispersion regime (λ>1.3μm), the opposite happens.

The GVD is different for different guided modes in MMFs. This is depicted in Fig. 3(h), where the value of β2(p)(ω) at 1.45 μm is shown for both SIMF and GIMF. The GIMF has smaller GVD values and exhibits narrower gaps across its mode groups compared to the SIMF, as further illustrated by the dispersion lengths plotted in Fig. 4(c). Consequently, GIMFs are expected to produce much less temporal broadening compared to SIMFs, suggesting that a smaller input power will be necessary to form multimode solitons in GIMFs.

By comparing the relevant characteristic lengths, we can find interesting conditions for the formation of solitons. For simplicity, here we focus on a fundamental soliton involving only one mode of an MMF. Similar to the case of singlemode fibers, a soliton involving only the mode p will have a hyperbolic secant profile: Ap(t)=P0(p)sech(t/T0),(17)where P0(p) and T0 are the peak power and duration of the soliton, respectively. From this, energy of this soliton is found to be E(p)=∫−∞∞|Ap(t)|2dt=2T0P0(p).(18)

By comparing the nonlinear length [Eq. (6)] with the dispersion length [Eq. (16)], we obtain the soliton order for mode p [56]: N(p)=LN(p)LD(p)=2πn2P0(p)T02λAeff(p)|β2(p)|.(19)

For N(p)<1, dispersive effects dominate over the nonlinear effects, preventing the formation of a soliton. When N(p)>1, higher-order solitons may form, exhibiting a periodic evolution. A fundamental soliton is formed when N(p)=1. Given that both the effective mode area (Aeff(p)) and dispersion (β2(p)) are different for different modes for an MMF, a soliton’s peak power will also depend on p. By setting N(p)=1, we find the following relations for the peak power P0(p) and the energy of a soliton of duration T0: P0(p)=λAeff(p)|β2(p)|2πn2T02,(20)E(p)=λAeff(p)|β2(p)|πn2T0.(21)

As the duration T0 varies, P0(p) and E(p) change as shown in Fig. 6(b) for the two types of fibers for the formation of solitons at λ=1550nm. For nanosecond pulses, the peak power needed is in the milliwatt range (or picojoule energy). For picosecond pulses, it reaches kilowatt levels (or nanojoule energy). For femtosecond pulses, the soliton needs peak powers in the megawatt range and energies reaching microjoule levels.

The dependence of the peak power on the mode number is shown in Fig. 6(a) for T0=105fs. The evolution of the required peak power closely matches the effective mode area in Fig. 3(e). Modes with larger effective mode areas require higher peak power (or energy) to form a soliton. The GIMF requires 4.6 kW peak power for mode 1 and 12 kW for mode 15, while the SIMF needs higher powers of around 22 kW for each mode.

The formation of an MMS inside an MMF requires compensation of temporal walk-off associated with different mode pairs. By equating the walk-off length LW(p,q) from Eq. (15) with the dispersion length in Eq. (16) and the nonlinear length in Eq. (6), we can determine the temporal duration of an MMS involving any two specific modes: Tw(p,q)=|β2(p)||β1(p)−β1(q)|.(22)

Based on the dispersion parameters of the two fibers in Fig. 3, Tw(1,6) of the first two radial modes is found to be 115 fs for a GIMF and 12 fs for a SIMF. Therefore, based on the walk-off compensation condition for forming an MMS, the duration and peak power for GIMF should be around 115 fs and 3.8 kW (see Section 2.B.1 where we provide numerical examples of MMS formation). Since β1(p) and β2(p) vary with wavelength [see Eqs. (7) and (8)], the values of Tw(p,q) also change. As a result, the “walk-off” soliton duration will also depend on wavelength (see Fig. 32 in Ref. [37]). For a soliton with T0>Tw(1,6), walk-off compensation becomes more critical, but it may be compensated for higher-order MMS with larger peak powers. For example, observing MMS in a SIMF requires a peak power in the MW range (or energy ∼100nJ) [42].

In the case of a soliton propagating in the fundamental mode, by substituting the mode effective area Aeff(p) with the beam area of the fundamental mode πR02 (R0 is the waist of the fundamental mode), one may rewrite the relation Eq. (21) between the soliton energy and duration as follows: M1=E(p)T0=λR02|β2(1)|n2.(23)

For MMSs, this relation can be approximated by replacing the mode energy E(p) with the total energy E, the mode effective area Aeff(p) with the beam area πRg2 (Rg is the beam waist) in an MMF, and the second-order dispersion β2(p) with that of the fundamental mode β2=β2(1) (supposing that the difference between chromatic dispersion of different modes can be neglected). Thus, the condition for the existence of an MMS can be expressed in terms of the beam size, pulse energy, and duration, as follows [62]: M≡ET0=λRg2|β2|n2>M1.(24)

As we shall discuss in Section 9, the experimental results have indeed demonstrated the validity of such a relationship (see Fig. 31 in Section 9).

Building on these foundational concepts, a thorough comparison of SIMF and GIMF can yield profound insights, facilitating a deeper understanding of MMS dynamics in multimode fibers.

An important consideration is that even when two pulses have identical durations and energies (or peak powers), their nonlinear lengths in the same MMF can differ vastly if their modal compositions are not the same. This suggests that the mode distribution of a pulse plays a crucial role in determining the nonlinear dynamics inside MMFs. For the same power and wavelength, the nonlinear length [Eq. (6)] of any mode is proportional to its MEA [Eq. (3)]. This feature indicates that modes with a smaller MEA have larger nonlinear coefficients [Eq. (5)] and they can excite a larger nonlinearity for the same given field power.

The MEA Aeff(p) in a GIMF increases with the mode order, while the MEA in a SIMF decreases, as seen in Fig. 3(e). This suggests that, for the same peak power, nonlinear effects may preferentially affect LOMs in GIMFs and HOMs in SIMFs, respectively. As an example, LOMs in a GIMF may initially compensate for their dispersion-induced pulse broadening via the Kerr effect, leading to the formation of an MMS composed by these modes [38].

Conversely, SIMFs are more likely to favor the formation of an MMS composed of HOMs, because of their smaller MEAs [see Fig. 3(e)] and higher MNCs. This is the reason that HOM Raman solitons are preferentially observed in SIMFs [42], while the Raman-beam cleanup is typically observed in GIMFs, rather than SIMFs [63]. With this insight, we can also understand Ref. [43], where a lower input-energy threshold was observed for generating a Raman MMS when the input beam was offset from the core’s center to enhance the excitation of HOMs of the fiber, compared to the case of on-center beam’s coupling (see Fig. 36 in Section 9.F).

In general, linear LOMs propagate faster than HOMs at the same wavelength, due to their larger group indices [see Fig. 5(b)]. In GIMFs, the Raman effect can initially cause a stronger red-shift for LOM components in individual pulses, because of their relatively higher nonlinear coefficients compared with HOMs. This red-shift decreases the speed of pulses in the LOMs [see the explanation of β1(p)(λ) in Section 2.C], thereby minimizing the group-velocity mismatch between LOMs and HOMs and facilitating the formation of an MMS.

Furthermore, for an MMS in GIMFs, the pulse in the fundamental mode is red-shifted the most, and HOMs with shorter wavelengths serve as pumps, thus facilitating a substantial energy transfer from HOMs to the fundamental mode via the Raman effect. This explains why a significant Raman beam clean-up is observed, as shown in Ref. [38]. This phenomenon is harder to occur in SIMFs because their fundamental mode does not have the largest nonlinear parameter.

The mode size increases in all MMFs as wavelength increases, which in turn decreases the MNCs. This feature potentially contributes to the observed decrease in the rate of red shift (SSFS) of a Raman MMS as it propagates over longer distances. Furthermore, at higher wavelengths the number of modes within an MMF diminishes, especially affecting HOMs, which lose their confinement as wavelength increases.

The nonlinear coupling among the modes of an MMF depends on the overlap of their mode profiles. Generally speaking, nonlinear coupling between two radial modes (l=0) tends to be stronger than that between a radial and a non-radial mode. Furthermore, within the same degenerate group, nonlinear coupling is stronger when compared with the coupling between two non-degenerate modes.

During mode competition within an MMF, modes with a larger MNC are more likely to dominate the nonlinear dynamics. This may partially explain the occurrence of self-cleaning observed in GIMFs [49]. It can also explain the observation of speckled patterns in SIMFs, which is indicative of the dominance of HOMs in SIMFs.

Building upon the insights gained from the preceding section, we discuss a widely used model for describing the evolution of optical pulses inside MMFs. As usual, the optical field of each pulse is written as the product of a rapidly evolving optical carrier exp[i(β0(0)z−ω0t)] and a slowly varying envelope E(x,y,z,t), where ω0 denotes the frequency at which the pulse’s spectrum peaks. The envelope E(x,y,z,t) is normalized in a way that |E(x,y,z,t)|2 is expressed in units of W/m2. This envelope can be expanded in terms of the eigenmodes of an MMF as follows: E(x,y,z,t)=∑p=1NFp(x,y)Ap(z,t),(25)where Fp(x,y) are the previously introduced transverse mode profiles and the mode amplitudes Ap have units of W.

The evolution of each modal component Ap along the propagation direction z is governed by the generalized multimode nonlinear Schrödinger equations (GMMNLSE) [64,65] ∂Ap(z,t)∂z=D{Ap(z,t)}+N{Ap(z,t)},(26)where the dispersion operator is given by D=i(β0(p)−β0(1))−(β1(p)−β1(1))∂∂t+i∑q≥24βq(p)q!(i∂∂t)q,(27)and the dispersion coefficients represent a Taylor-series expansion of β(p)(ω) in Eq. (7) evaluated at the frequency ω0: $$\begin{gathered} β(p)(ω)=β0(p)+β1(p)(ω−ω0)+12β2(p)(ω−ω0)2 \\ +∑n≥3∞1n!βn(p)(ω−ω0)n. \end{gathered}$$(28)

It is assumed that the pulse evolves in a reference frame moving at the velocity of the fundamental mode 1/β1(1). All phase and group velocities are referenced to the corresponding values for the fundamental mode, i.e., the first two terms in Eq. (27) are relative to the fundamental mode.

The modes are coupled by the nonlinear effects such as SPM, XPM, and FWM, whose impact in Eq. (26) is included by the following nonlinear terms: $$\begin{gathered} NA=in2ω0c(1+fSiω0∂∂t)∑l,m,nN[(1−fR)SplmnKAlAmAn* \\ +fRSplmnRAl∫−∞thR(τ)Am(z,t−τ)An*(z,t−τ)dτ], \end{gathered}$$(29)where SplmnK and SplmnR are the nonlinear mode-coupling coefficients for the Kerr and Raman effects, respectively [64]. The values of SplmnK and SplmnR are in general different when the modes’ polarization is taken into account, whereas they become identical within the scalar description of the fiber modes that we adopt in this section [66]. Therefore, we can drop the superscripts K and R and express the nonlinear coefficients as Splmn=∬dxdyFpFlFmFn∬dxdyFp2∬dxdyFl2∬dxdyFm2∬dxdyFn2,(30)where we made the further assumption of real-valued transverse mode profiles. Equation (29) incorporates the effects of Kerr nonlinearity, Raman scattering (associated with the coefficient fR), and self-steepening (related to the coefficient fS). The nonlinear index of silica is n2=2.7×10−20m2/W, and its Raman fraction is fR=0.18. The Raman response function, which is crucial for capturing the nonlinear response of an MMF to ultrashort pulses, can be approximated by hR(t)=(τ1−2+τ2−2)τ1exp(−t/τ2)sin(t/τ1),(31)where τ1=12.2fs and τ2=32fs are related to the frequency and the damping time of the Raman response in silica [56].

This subsection provides concise examples illustrating various phenomena, such as the linear multimode walk off, formation of multimode solitons, their fission and creation of Raman solitons, and collision of two such solitons. These examples not only serve to clarify the basic concepts but also provide a basis for both analytical and experimental exploration.

We assume that the input conditions for all modes are the same in the form of Gaussian pulse, centered at t0 with a duration characterized by its full width at half maximum (FWHM) TFWHM=1.665T0. More specifically, Ap(z=0,t)=CpET0πe−(t−t0)22T02,(32)where Cp represents the initial modal coupling coefficient, and E is the input energy. They are found by projecting an input Gaussian beam Ein(x,y) onto the eigenmode basis, expressed mathematically as Cp=∬−∞∞Ein(x,y)·Fp(x,y)dxdy.(33)

Input conditions are varied by adjusting the two width ratios for the input beam. They are defined as $$\begin{gathered} rw=w/wmode,p=1, \\ rs=x0/wmode,p=1, \end{gathered}$$(34)where w is beam’s width, x0 is its lateral offset, and wmode,p=1 is the width of the fundamental mode. The use of these dimensionless parameters allows for describing diverse input conditions for a given input pulse with parameters such as input energy E, temporal duration T0, and the parameter values rw and rs.

When an input beam is aligned with the center of the fiber’s core, only a few radial modes are excited (specifically modes p=1,6,15,…, as depicted in Fig. 3). For simplicity, we only take into account these three dominant radial modes corresponding to p=1,6,15 in the following examples.

In the linear regime, modal dispersion leads to temporal walk-off of pulses carried by different modes in a multimode GIMF; see Figs. 7(a) and 7(b). In this case, an input beam with rw=3 can be decomposed into the first three radial modes (p=1, 6, and 15). The input pulse has a Gaussian shape with a duration of T0=110fs (FWHM: 184 fs) and has an energy of E=0.01nJ. As the peak power of the input pulse is merely 27.5 W, nonlinear effects remain insignificant due to LN(1)≈52.2m, but dispersion effects play a pivotal role. This can be seen by noting that the dispersion length is LD(1)≈0.72m. The walk-off lengths LW are 0.82 m and 0.40 m for the mode pairs (1.6) and (1,15), which correspond to DMDs of 0.224 ps/m and 0.464 ps/m, respectively. As seen in Fig. 7(a), after propagating over a distance of 27.7LD=20m, the multimode pulse undergoes significant temporal broadening. The peak delay between mode 1 and mode 6 is approximately 4.8 ps and it increases to 9.3 ps for modes 1 and 15. Given that dispersion solely alters the spectral phase, the spectral power distribution remains unaffected, as illustrated in Fig. 7(b).

However, when the input energy is increased to 0.7 nJ, corresponding to a peak power P0=3.58kW, the peak powers of the three excited modes become 1.76 kW for mode 1, 1.10 kW for mode 6, and 0.72 kW for mode 15, and result in different nonlinear lengths for these modes: LN=0.82m for mode 1, 2.55 m for mode 6, and 5.62 m for mode 15. Although N(1)=0.94<1 without considering higher-order mode powers, a qualitatively distinct phenomenon is observed, as shown in Figs. 7(c) and (d). Mode 1 experiences significantly less temporal broadening, facilitating the trapping of a substantial portion of the energy from mode 6. This interaction leads to the formation of an MMS that evolves as a single entity. The group delay of the MMS in Fig. 7(c) falls in between that of modes 1 and 6 in Fig. 7(a). This feature shows that an averaging of the linear group velocities of the constituent modes occurs when an MMS forms. In order to sustain this balanced group velocity, the Kerr nonlinearity produces a red-shift for mode 1 (which decreases its speed) and a blue-shift for mode 6 (which increases its speed) [see Fig. 7(d)]. Note also that the peak power of mode 15 is not large enough to be trapped into the MMS. As a result, the pulse in this mode spreads its energy into a dispersive wave, without much change in its spectral characteristics, as seen in Fig. 7(d).

By varying the input pulse energy, we can determine the threshold for MMS formation. In Figs. 8(a)–8(c), we plot the temporal power distribution (normalized to its maximum value) for modes 1, 6, and 15 at the fiber’s output as a function of input pulse energy. In the quasi-linear regime (0<E<0.5nJ), three modes exhibit significant walk-off and pulse broadening. However, owing to its large nonlinear coefficient γ(1) and its high peak power, mode 1 begins to exhibit less broadening as the input energy increases. To illustrate this feature quantitatively, we plot in Figs. 8(d)–8(f) the pulse peak temporal position, the peak power, and the pulse duration. The temporal position of peak power for mode 1 exhibits a large delay when an MMS forms by trapping mode 6 [see Fig. 8(d)], while the power of all modes keeps increasing [see Fig. 8(e)].

However, only the duration of mode 1 is significantly reduced [see Fig. 8(f)]. At an energy of around 0.5 nJ, the peak position of mode 6 aligns with mode 1 [see Fig. 8(d)], which can be considered as the threshold for MMS generation. For mode 15, even though a bit of energy is trapped, most of its energy flows away as a dispersive wave [see Fig. 8(c)], until 0.8 nJ, which is also observed as the turning point of peak power changes in Fig. 8(e). As the MMS forms, the durations of pulses in all three modes become very close to each other, as shown in Fig. 8(f).

Next, let us consider what is the influence of the input conditions (specifically, the input modal composition) on the formation and propagation characteristics of an MMS, for a given input energy. In Figs. 9(a)–9(c) we plot the output temporal profiles of power (normalized to their maximum values) in modes 1, 6, and 15, respectively, versus input beam waist width rw. Here the input energy is fixed at E=1nJ. The corresponding pulse peak temporal position, the pulse peak power, and the pulse duration are depicted in Figs. 9(d)–9(f), respectively. At this energy level, the input pulse energy is effectively redistributed among modes 1 and 6 in order to efficiently form an MMS (see Fig. 8), whereas a portion of input energy coupled to mode 15 only contributes to dispersive waves.

When rw=1, only mode 1 is excited, since the input beams have the same size as mode 1. A singlemode soliton within mode 1 is formed, as indicated by the zero peak power for modes 6 and 15. As rw increases, the group velocity of the MMS decreases, which is reflected in the larger group delay at the output, as shown in Fig. 9(d). Moreover, with an increase in rw, the high-order modal content of the MMS in modes 6 and 15 also increases, as can be seen in Fig. 9(e). The duration of the MMS expands from 0.1 ps to 0.17 ps, wherein the duration of the pulse in the fundamental mode is slightly shorter than that of pulses in the higher-order modes, as demonstrated by Fig. 9(f).

In the previous examples, the Raman effect was not considered (fR=0) in the GMMNLSE model. Now, we aim to illustrate how the Raman effect affects MMS. To achieve this, we introduce the Raman effect by setting the delayed nonlinear response coefficient to fR=0.18 in Eq. (29). Once the input energy exceeds 1 nJ, the Raman effect becomes notably influential. Here, we provide two examples to showcase its impact on MMS.

In Fig. 10 we illustrate the evolution of the MMS for an input pulse with E=3nJ and a duration of T0=110fs (FWHM: 184 fs), propagating through a 2 m GIMF. At this energy level, nonlinear effects intensify, enhancing mode coupling and inducing significant temporal and spectral changes. The LOMs undergo a red-shift, whereas the HOMs experience a blue-shift at z=0.4m [see Fig. 10(b)], in order to maintain a GVD-induced balance of their different group velocities. As a result, a Raman MMS forms, with a relatively narrower temporal duration and a noticeably red-shifted spectrum, owing to SSFS. At z=2m, the central spectrum of the Raman MMS shifts to 1.52 μm. This is accompanied by a decrease in the soliton group velocity. When comparing with the MMS group delay in Fig. 7(c), where a 20 m propagation leads to approximately 1.5 ns delay, the MMS in Fig. 10(a) only propagates over 2 m.

When increasing the input energy up to 7 nJ, as depicted in Figs. 10(c) and 10(d), there is sufficient nonlinearity to reorganize the field around t=0 into a second MMS. This suggests that the initial high-order soliton undergoes fission, resulting in two different MMSs. The first soliton, which carries most of the energy, exhibits a greater delay (extending beyond the plotting window at 1.6 m and 2.0 m) and a larger SSFS. The second MMS carries less energy; however, this is sufficient to form another Raman MMS, featuring a reduced delay and SSFS. Further energy increments can lead the second MMS to gather more energy, resulting in larger group delays (slower group velocities) and heightened SSFS, as detailed in the subsequent subsection.

In previous subsections, we observed different dynamics of multimode solitons in MMFs, including modal walk-off, formation of an MMS, and its fission. We also discussed the dependence of the MMS properties on physical parameters such as the input pulse’s duration, its group velocity, and SSFS upon varying its energy or modal content. Here, we discuss the possibility of controlling the SSFS of an MMS by adjusting its modal content. Under specific conditions, two Raman MMSs with different group velocities may undergo inelastic collisions, resulting in an energy exchange between them.

As an example, we cite the result in Ref. [39]. In this study, the input beam led to the excitation of the first 15 modes, and the input pulse’s duration was 70 fs with a center wavelength of 1.4 μm. Figure 11(a) shows simulated output spectra for a 2 m long GIMF as a function of input pulse energy using rw=2 and rs=0, corresponding to the excitation of only radial modes. The SSFS of the first MMS (S1) at the output is visible for input energies exceeding 1 nJ, while the SSFS of the second MMS (S2) becomes pronounced beyond 5 nJ. The SSFS of both Raman MMSs increases with growing input energy.

When the input beam is slightly shifted from the center of the fiber axis, additional non-radial HOMs are introduced (rw=2.8, rs=0.5). In this case, a similar spectral evolution occurs as seen in Fig. 11(b). However, at 15 nJ, S1 undergoes less SSFS-induced red-shift and temporal delay when compared with S2, owing to its reduced fundamental mode content. As the input beam undergoes further energy increase and axial position shift (rw=3.2, rs=1.7), the special spectral evolution of S1 seen in Fig. 11(c) becomes evident. The corresponding experimental demonstration is shown in Fig. 40 of Section 9.H. As seen in Fig. 11(c), a surprising spectral feature occurs for input energies between 22.5 nJ and 26 nJ: the soliton spectra exhibit large separations before merging again.

To understand these dynamics, we show the temporal and spectral evolutions of the MMS in Figs. 11(d) and 11(e): here the input energy is 22.7 nJ. Initially, the input pulse splits into two solitons (S1 and S2), which collide and then separate due to their inelastic energy exchange. Despite S2 having a smaller peak power when compared with S1 before their collision [see Fig. 11(f)], S2 gradually exhibits a larger delay (and larger SSFS), eventually leading to its collision with S1.

The analysis of the modal contents reveals that S2 has a relatively larger contribution from the fundamental mode (S2M1) when compared with S1 (S1M1), as depicted in Fig. 11(g). This occurs even though the sum of other mode contents in S1 (S1M2-15) is larger than that of S2 (S2M2-15), as shown in Fig. 11(h). Since the fundamental mode dominates the nonlinear strength, S2 experiences a larger SSFS and group delay, resulting in its collision with S1. As a result of their inelastic collision, the two MMSs exchange energy and separate, because S1 gains more energy in the fundamental mode, thus exhibiting a larger group delay. In Section 9, we shall describe observations giving experimental evidence of the spectral evolution and associated dynamics of collisions between MMSs in multimode GRIN fibers [39].

Furthermore, MMS interactions have been theoretically studied in several papers. In Ref. [67], the interaction between two MMSs in a GIMF, each carrying only the degenerate modes LG11a and LG11b with consideration of their polarization states, was studied. In analogy with the case of single-mode fibers, these MMSs exhibit an attraction or a repulsion behavior, on the basis of their relative phase, despite belonging to different fiber modes. However, unlike the single-mode scenario, each soliton transfers some power to the other mode via intermodal four-wave mixing. Despite this intermodal power exchange, both MMSs maintain their bimodal nature, and interact as if they were traveling in a single-mode fiber.

In Ref. [68], the interaction between fundamental solitons in single-mode fibers under the influence of Raman scattering was extensively investigated. The study reveals a net energy transfer from the leading pulse to the trailing pulse, akin to the behavior that we have described for the multimode scenarios. Furthermore, the impact of different soliton amplitudes or phases on the fraction of transferred energy was analyzed.

In this section, we first introduced the eigenmodes and eigenvalues of MMFs, which are important for understanding the nonlinear as well as linear properties of an MMF. By comparing two specific MMFs, a SIMF and a GIMF, we found that the GIMF has dominant nonlinear coefficients for LOMs, thus favoring MMSs composed of low-order modes. Moreover, a GIMF has much reduced modal dispersion, which greatly facilitates the MMS formation. In contrast, a SIMF favors the generation of HOM solitons.

Second, by analyzing the modal eigenvalues, we obtain a piece of key information about the propagation properties of MMFs. For instance, equally-spaced propagation constants in GIMFs enhance mode beating, giving rise to the self-imaging effect. Thus, analyzing β1 for different modes and wavelengths helps us to better understand a key property of the MMS in GIMFs: LOMs propagate faster than HOMs in linear regimes. Therefore, the overall group velocity of MMSs is influenced by their modal content. Building on these fundamental notions, we introduced the GMMNLSE, which is a relevant model for describing ultrashort pulse propagation in a broad range of MMFs with identifiable eigenmodes.

We elucidated a series of phenomena typical of GIMF: MMS formation via compensation of modal walk-off, MMS fission, mode-based SSFS control, and MMS collisions. Specifically, we analyzed how the Raman effect interacts with MMS: it increases the soliton’s wavelength, and LOMs exhibit larger SSFS in GIMFs. This has two significant impacts: first, it reduces the velocity of LOMs, thus diminishing their group velocity mismatch with HOMs and favoring their nonlinear trapping. Second, LOMs with increased red shifts are fed energy by the HOMs with a reduced red-shift, resulting in significant energy transfer from HOMs to LOMs, and eventually leading to Raman beam cleanup into the fundamental mode: corresponding experiments are discussed in Section 9.

Lastly, we discussed the possibility of SSFS manipulation by changing the input modal composition of an MMS. Under specific conditions, two MMSs resulting from the fission of an input pulse and propagating with different modal contents may collide at some point along the fiber, owing to their different group velocities. This section lays the groundwork for deeper explorations into MMS, by either analytical methods or experiments, as discussed in subsequent sections.

In this section, we will review the theoretical framework, leading to analytical results, for describing the dynamics of MMSs. Most of these results rely on the variational formulation of Eq. (26) and the derivation of effective theories for describing an MMS by means of a reduced number of degrees of freedom. In this context, most of the literature focuses on two-mode (or bimodal) soliton formation [69–72]. The multimode case, whose treatment is much more complex, has received much less attention [69]. In Section 4.A, we introduce first the variational formulation for the two-mode soliton case. The Kantarovitch, or average variational approach, which is our main method in this section, is presented in Section 4.B. The effective variational theory obtained with this method for the case of a two-mode bright soliton is presented in Section 4.C. Later, in Section 4.D we will present the virial, or moment method, reviewing the results presented in Ref. [72]. Finally, in Section 4.E we present a generalization of this analysis to a multimode context, following Ref. [69].

When just considering two modes, and neglecting high-order effects (e.g., high-order dispersion, Raman effect), Eq. (26) can be written as $$\begin{gathered} ∂zA1=iΔβ0(1)A1−Δβ1(1)∂tA1+iβ2(1)2∂t2A1 \\ −i(R11|A1|2+2R12|A2|2)A1, \end{gathered}$$(35a)$$\begin{gathered} ∂zA2=iΔβ0(2)A2−Δβ1(2)∂tA2+iβ2(2)2∂t2A2 \\ −i(R22|A2|2+2R21|A1|2)A2, \end{gathered}$$(35b)where we have considered the tensor Rplmn=Splmnn2ω0/c symmetries [35] Rmmnn=Rmnmn=Rmn.(36)

After a straightforward normalization (see, for example, Ref. [73]), Eq. (35) reduces to i(∂z+δ∂t)u+νuu+βu2∂t2u+(|u|2+C|v|2)u=0,(37a)i(∂z−δ∂t)v+νvv+βv2∂t2v+(C|u|2+|v|2)v=0,(37b)with βu,v=1 for anomalous and βu,v=−1 for normal GVD. Note that, under the phase transformation u→eiνuzu(z,t),v→eiνvzv(z,t),(38)Eq. (37) reduces to i(∂z+δ∂t)u+βu2∂t2u+(|u|2+C|v|2)u=0,(39a)i(∂z−δ∂t)v+βv2∂t2v+(C|u|2+|v|2)v=0,(39b)which is the starting point for most of the analytical research that has been performed in this system (see, e.g., Refs. [70,71]). A key parameter here is the cross-phase modulation (XPM) coefficient C, which couples fields u and v. With C=1, Eq. (39) is known as the Manakov equations [74]. Alternatively, one can get rid of the walk-off terms by performing the transformation u(z,t)=U(z,t)eiδ2z/2−iδt,V(z,t)=v(z,t)eiδ2z/2+iδt,(40)which yields the system i∂zU+βu2∂t2U+(|U|2+C|V|2)U=0,(41a)i∂zV+βv2∂t2V+(C|U|2+|V|2)V=0.(41b)

In this context, the effect of the walk-off can be reintroduced as an initial splitting velocity between the modal components of a pulse as it is injected into the fiber. Thus, the effect of temporal walk-off owing to modal dispersion can be transferred to an initial relative frequency shift among the two input pulse components [75,76].

The set of Eq. (37) can be derived from the Lagrangian density L≡Lu+Lv+Lint,(42)with $$\begin{gathered} Lu≡i(u*∂zu−u∂zu*)+iδ(u*∂tu−u∂tu*)−βu|∂tu|2+|u|4, \\ Lv≡i(v*∂zv−v∂zv*)−iδ(v*∂tv−v∂tv*)−βv|∂tv|2+|v|4, \end{gathered}$$and the interaction contribution Lint≡2C|u|2|v|2.

Then Eq. (39) can be retrieved through the Euler-Lagrange equations ∂∂z(∂L∂uz*)+∂∂t(∂L∂ut*)−∂L∂u*=0,(43a)∂∂z(∂L∂vz*)+∂∂t(∂L∂vt*)−∂L∂v*=0.(43b)

Similarly, the application of the transformation Eq. (40) to the Lagrangian Eq. (42) leads to the new Lagrangian density L≡LU+LV+Lint,(44)with LU≡i(U*∂zU−U∂zU*)−βu|∂tU|2+|U|4,(45)LV≡i(V*∂zV−V∂zV*)−βv|∂tV|2+|V|4,(46)Lint=2C|U|2|V|2,(47)from where Eq. (41) can be also derived.

By defining the generalized momenta Pu≡∂L∂Uz∗=iU/2,Pv≡∂L∂Vz∗=iV/2,(48)this system can be alternatively described by Hamiltonian density obtained through the Legendre transform H=PuUz∗+Pu∗Uz+PvVz∗+Pv∗Vz−L,(49)which leads to H=HU+HV+Hint,(50)with HU=βu|∂tU|2−|U|4,HV=βv|∂tV|2−|V|4,(51)and Hint=−2C|U|2|V|2.(52)

From the momentum-energy tensor defined through our Lagrangian [77], the following conservation laws can be deduced. •Conservation of power. The powers carried by each mode are defined as Pu≡∫−∞∞|u|2dt,Pv≡∫−∞∞|v|2dt.(53)These quantities satisfy ∂zPu=0,∂zPv=0,(54)which means that they are conserved individually along propagation in z. This is a direct consequence of the fact that all terms in Eq. (42) are real [72].•Conservation of the total momentum. The momentum for each field u, v is defined as Ju≡i∫−∞∞(uut*−u*ut)dt,Jv≡i∫−∞∞(vvt*−v*vt)dt.(55)The total momentum J=Ju+Jv is conserved along propagation, i.e., ∂zJ=∂z(Ju+Jv)=0.(56)•Conservation of the total energy. The last quantity to be conserved is the total energy or Hamiltonian function H=Hu+Hv+Hint=∫−∞∞(Hu+Hv+Hint)dt.(57)The conservation law in this case reads ∂zH=∂z(Hu+Hv+Hint)=0.(58)

In this section, we introduce the averaged variational method, which is widely used in order to compute soliton solutions in an optical context [15]. Originally presented in an optical context by Anderson et al. [78], this method is also known as the Rayleigh-Ritz or Kantarovitch optimization approach and generalizes its stationary version: the Ritz approach [78–80]. This method allows us to compute approximate analytical solutions of a given nonlinear partial differential equation by applying the principle of least action to a parameter-dependent solution ansatz, based on either the Lagrangian or the Hamiltonian formalism. In what follows, we describe the main steps in either the Lagrangian or the Hamiltonian formalism.

The method consists of the following four main steps. (1)First of all, we need to define an approximate ansatz solution, or trial function, that captures the main features and shape of the state that we want to compute. This ansatz is a function of the form u=u[t;q(z)], which depends on z through a number of parameters: q(z)={q1(z),…,qn(z)},which are the generalized coordinates of the system.(2)Given a Lagrangian density associated with a certain set of partial differential equations [see, e.g., Eq. (42)], the next step in this procedure is to compute the effective reduced Lagrangian function of the system, which, in the context of temporal MMSs studied here, is defined as L[q(z)]≡∫−∞∞L[u,u*,…;q(z)]dt.(59)(3)After that, we obtain the dynamical system for the generalized coordinates by computing the Euler-Lagrange equations ddz(∂L∂(dzqm))−∂L∂qm=0(60)for each parameter qm with m=1,…,n, where we have defined dzqm≡dqm/dz.(4)Finally, in the last step, we study the dynamics of the reduced system Eq. (60), and rebuild the behavior of the original system by using the initial ansatz.

Equivalently, one may consider the Hamiltonian formalism for obtaining an effective reduced system. The process is equivalent to that followed in the Lagrangian case, but now we use the Hamiltonian H[q(z),p(z)]≡∫−∞∞H[u,u*,…;q(z)]dt,(61)where p(z)={p1(z),…,pn(z)} are the generalized momenta defined as pm=∂L∂(dzqm).

Then, the effective dynamics of the system are captured by the Hamiltonian equations of motion dqmdz=∂H∂pm,dpmdz=−∂H∂qm(62)for each m=1,…,n.

Once that variational machinery has been introduced, we will try to answer the following question: which kind of initial condition for Eq. (39) leads to the formation of a two-mode soliton?

Most of the works dealing with two-mode soliton formation focus on a regime where each mode in the fiber experiences anomalous GVD [69–71,75]. In comparison with this case, the mixed scenario has been much less studied [81].

In this section, we will mainly focus on the interaction of bright solitons, where both modal components are propagating in the anomalous dispersion regime. In what follows, we take βu=βv=1. In this context, two main behaviors are found: either (i) the two solitons split and escape or (ii) they are trapped, forming a bound state, i.e., a two-mode soliton. These solitons, however, are not static upon propagation, but they undergo oscillations in their separation and peak powers [73]. An example of this regime is depicted in Fig. 12.

The threshold between these two regimes deeply depends on two main factors: the strength of temporal walk-off, or the initial relative frequency shift of the pulses, and the strength C of the XPM. One of the first researchers who investigated these phenomena was Menyuk in 1987. He found, by performing numerical simulations, that the splitting could be arrested through XPM, and that the threshold for mutually trapped solitons increases with increasing walk-off [73]. After that, many theoretical works focused on understanding this phenomenon and on the analytical determination of that threshold [69–72].

Instead of reviewing these results in chronological order, we have decided to organize them in terms of the variational formalism that was applied.

Let us first review some results based on a Lagrangian approach. One of the first works to tackle this problem was that of Ueda and Kath [70]. In the following, we will use these results as the main guiding example for didactic purposes.

Their approach considered the Lagrangian Eq. (44), which is more adequate for investigating soliton formation, in detriment of the explicit appearance of the linear walk-off δ in the results.

The general soliton ansatz they considered was u=ηusech(t−τuwu)eiΦu,v=ηvsech(t−τvwv)eiΦv,(63)with the phases $$\begin{gathered} Φu≡Vu(t−τu)+θu2wu(t−τu)2+ϕu, \\ Φv≡Vv(t−τv)+θv2wv(t−τv)2+ϕv, \end{gathered}$$where ηu,v, wu,v, and τu,v are the amplitude, width, and central position of the fields u and v, respectively. These parameters correspond to the three lowest-order moments of the soliton envelopes. Vu,v, θu,v, and ϕu,v, represent the three lowest-order moments of the phase Φu,v. These are the nonlinear velocity of the soliton central position, the chirp of the soliton (related to the change of the soliton temporal width), and the linear phase, respectively.

By computing the Euler-Lagrange Eq. (60) with this ansatz, the authors obtained an eight-dimensional reduced dynamical system for describing the evolution of the corresponding ansatz parameters. This effective theory allowed the study of both symmetric and asymmetric soliton configurations. However, due to the complexity of the former scenario, the authors focused on the symmetric case, which is the one that we review in the following paragraphs.

Let us consider two identical solitons that are symmetrically located in time on opposite points around their mean position, and that propagate in opposite directions (in a reference frame propagating with the average group velocity). This implies: wu=wv=w, ηu=ηv=η, θu=θv=θ, Vu=−Vv=−V, and τu=−τv=−τ. With these considerations, the effective Lagrangian L=∫−∞∞L(u,u*,ν,ν*)dt(64)yields the Euler-Lagrange equations dτdz=V,(65a)dVdz=CPw2∂ΔF(Δ),(65b)dwdz=θ,(65c)dθdz=4π2w2{w−1−P2−32CP∂Δ[ΔF(Δ)]},(65d)where F(Δ)≡coshΔsinh3Δ(Δ−tanhΔ),Δ≡2τw,(66)and P=Pu=Pv=2η2w is the conserved mode power [see Eq. (53)]. The fixed points of this system correspond to a shape-preserving two-mode soliton, which is a center of the dynamics; hence it is neutrally stable: any small fluctuation of the system will drive it away from this point [82]. Interestingly, around this point a continuous family of oscillatory states may appear, which correspond to dynamical two-mode solitonic states. We will come back to this in a few paragraphs.

A way to obtain more information about this bimodal soliton system is by studying the effective Hamiltonian H=Hu+Hv+Hint=∫−∞∞(HU+HV+Hint)dt,(67)associated with our ansatz. After some algebra, the Hamiltonian can be reorganized as H=HI+HK, where HI≡16(P2−1w)2+π224θ2(68a)is a measure of the internal energy in each soliton, which is taken up by pulse-width oscillations and HK≡12V2−CP2wF(Δ)(68b)represents the kinetic and potential energy due to the coupling between the u and ν solitons [70]. As we know from the conservation law Eq. (58), H=E=const, where E is the initial energy of the system, that is, E≡H(0).

Depending on the sign of E, two behaviors are possible. •E<0: bound state regime. The two pulses attract each other, and form a bound state, wherein the pulse central position τ and width w undergo periodic oscillations in z [70,73]. This bound state corresponds to a dynamic two-mode soliton.•E>0: escaping or splitting regime. When H=E>0, however, the two pulses have enough initial velocities, or frequency-shift (walk-off) to overcome the XPM-induced attraction, yielding two well-separated independent solitons.