Introduction

The past few years have witnessed a resurgence of interest in multimode waveguide structures, predominantly driven by the ever-increasing demand for higher information capacities1,2,3,4. This renaissance, in turn, incited a flurry of activities in the general area of nonlinear multimode fiber optics5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25. In this regard, the sheer complexity associated with the presence of hundreds or thousands of nonlinearly interacting modes that collectively act as a many-body system, has led to new opportunities in observing a multitude of novel optical effects that would have been otherwise impossible in single-mode settings12. These include for example, the formation of multimode solitons (MM-solitons)11, beam self-cleaning effects16,17,19, geometric spatiotemporal instabilities14, and efficient supercontinuum and second-harmonic generation17,20,21,25, to mention a few. The prospect for spatiotemporal laser mode-locking has also been successfully demonstrated in a recent study26.

One of the most prominent processes taking place in nonlinear optical waveguides is that associated with dispersive wave (DW) emission or Cherenkov radiation27,28,29,30,31,32. In general, these effects are induced by optical solitons generated in the anomalous dispersion regime. In single-mode fibers (SMFs), this Cherenkov emission is typically the outcome of a soliton fission event, caused by the inherent instability of higher-order solitons. This latter highly nonlinear effect significantly expands the spectrum, thus allowing phase-matching to occur between the soliton itself and the expected Cherenkov lines. In turn, this condition promotes the generation of a DW in the normal dispersion region of a SMF. Nonlinear interactions between DWs, solitons and other spectral features, are known to play a key role in shaping supercontinuum generation, especially in single-mode, dispersion-engineered photonic crystal fibers and integrated microstructures33,34,35.

As one could expect, this situation is considerably more involved in nonlinear multimode environments. In fact, in a recent experimental study12, carried out in a parabolic-index multimode fiber (MMF), a series of dispersive wave lines was unexpectedly observed around the visible wavelengths. In this case, this Cherenkov comb extended over hundreds of nanometers, as opposed to SMFs where DW phase-matching can only be observed within a narrow band. Naturally, the question arose as to what leads to the creation of such distinct DW lines. A possible explanation for this effect was provided in ref. 18, by interpreting these lines as Kelly sidebands—a characteristic sideband instability associated with periodically amplified solitons36. After all, in a parabolic fiber the soliton beam tends to undergo fast periodic compressions and expansions, thus mimicking nonlinear effects emanating from periodic amplification, even in a fully passive MMF. What makes this periodic self-imaging behavior possible is the equidistant distribution of the propagation eigenvalues, something that is very unique to parabolic (harmonic) potentials. While this theoretical model could elucidate some of the trends in the DW spectrum, it could not account for many of the observed spectral feature. Even more puzzling was the appearance of blurred DW bands in experiments performed in step-index MMFs, where no such periodic beam compression/expansion is possible. This directly suggests that the physics behind the MMF Cherenkov lines is considerably more complex than previously thought—implying that the actual mechanism behind this effect is still elusive.

In this article, we develop a theoretical model capable of predicting and explaining the Cherenkov spectral peaks generated by multimode soliton fission processes in nonlinear MMFs. A key notion in our model is the very formation of a multimode soliton, that forces all the modes involved to first coalesce in the temporal domain, and in doing so, individually shift their spectral content. The actual Cherenkov lines then ensue from a wideband multimode phase-matching with these MM-solitons after fission and, in principle, they can lie in the anomalous dispersive region, in contradistinction to SMFs. Our theoretical model is in close agreement with experimental observations in both parabolic and step-index MMFs. In addition to enabling promising applications in nonlinear MMFs, our results can shed light on these highly involved nonlinear many-mode effects.

Results and discussion

Theory

A central element in our work is the way a MM-soliton is established when excited in the anomalous dispersive region of a MMF. In general, a multimode soliton is a composite entity in which several modes bond and hence travel together at the same velocity V. To form a multimode soliton, the nonlinearity should first compensate for both intermodal dispersion (differential group delays) and chromatic dispersion. Even more importantly, the carrier frequencies ωμν associated with each mode must be nonlinearly shifted with respect to each other so as to lock all the modes at the same group speed V. For simplicity, we here ignore any polarization effects, in which case the nonlinear contribution to the refractive index can be written as \(n\left( {{\boldsymbol{r}},z,\omega } \right) = n_0\left( {{\boldsymbol{r}},\omega } \right) + n_2\left| {\hat E\left( {{\boldsymbol{r}},z,t} \right)} \right|^2\), where n0 (r, ω) represents the waveguide index profile and n2 the nonlinear Kerr coefficient of silica glass. Following the above arguments, we express the total electric field \(\hat E\left( {{\boldsymbol{r}},z,t} \right)\) in this system as a superposition of all the spatial modes involved, i.e.,

$$\hat E\left( {{\boldsymbol{r}},z,t} \right) = \mathop {\sum}\limits_{\mu ,\nu } {E_{\mu \nu }\left( {\boldsymbol{r}} \right)} \exp \left[ {i\beta _{\mu \nu }\left( {\omega _{\mu \nu }} \right)z - i\omega _{\mu \nu }t} \right]{\Phi}_{\mu \nu }\left( {z,t} \right),$$
(1)

Where Eμν (r) denotes the modal field profile of mode (μ, ν), and βμν, ωμν represent its corresponding propagation constant and angular frequency. Meanwhile, Φμν (z,t) stands for the slowly varying envelope associated with the (μ,ν) mode, where μ,ν {0, 1, 2, 3, …}. By following the analysis of Crosignani and Di Porto37,38, one then finds that, in the anomalous dispersion region, the envelopes evolve according to

$$\, i\left( {\frac{{\partial {\Phi}_{\mu \nu }}}{{\partial z}} + \frac{1}{{v_{\mu \nu }}}\frac{{\partial {\Phi}_{\mu \nu }}}{{\partial t}}} \right) + \frac{1}{2}\left| {\beta _{\mu \nu }^{\left( 2 \right)}} \right|\frac{{\partial ^2{\Phi}_{\mu \nu }}}{{\partial t^2}}\\ \, + \mathop {\sum}\limits_{{ {\mu {\prime} \ne \mu } \atop {\nu {\prime} \ne \nu } }} {\left[ {R_{\mu \nu ,\mu \nu }\left| {{\Phi}_{\mu \nu }} \right|^2 + 2R_{\mu {\prime}\nu {\prime},\mu \nu }\left| {{\Phi}_{\mu {\prime}\nu {\prime}}} \right|^2} \right]} {\Phi}_{\mu \nu } = 0,$$
(2)

where \(\beta _{\mu \nu }^{\left( 2 \right)} = \left( {d^2\beta _{\mu \nu }/d\omega ^2} \right)_{\omega = \omega _{\mu \nu }},v_{\mu \nu }^{ - 1} = \left( {d\beta _{\mu \nu }/d\omega } \right)_{\omega = \omega _{\mu \nu }}\) represent the second-order dispersion coefficient and modal-group velocity. The self-phase and cross-phase modulation coefficients are given by \(R_{\mu {\prime}\nu {\prime},\mu \nu } = \left( {\omega n_2/c} \right)\int\!\!\!\int E_{\mu {\prime}\nu {\prime}}^2\left( {\boldsymbol{r}} \right)E_{\mu \nu }^2\left( {\boldsymbol{r}} \right)d{\boldsymbol{r}}/\int\!\!\!\int E_{\mu \nu }^2\left( {\boldsymbol{r}} \right)d{\boldsymbol{r}}\)37,38. In developing Eq. (2), we intentionally omitted any four-wave mixing (FWM) terms. This FWM suppression can be justified37,38 given that the frequencies ωμν of the participating modes will be eventually quite different once a MM-soliton is formed. As previously noted, a MM-soliton can only form provided that all the constituent modes move at a common group velocity V, i.e., νμν = V. In this case, one can show that Eq. (2) can admit the following self-consistent multimode bright-soliton solution,

$${\Phi}_{\mu \nu }\left( {z,t} \right) = {\Phi}_0^{\mu \nu }\exp \left( {\frac{{i\left| {\beta _{\mu \nu }^{\left( 2 \right)}} \right|z}}{{2\tau ^2}}} \right){\mathrm{sech}}\left( {\frac{{t - z/V}}{\tau }} \right),$$
(3)

where \({\Phi}_0^{\mu \nu }\) is the peak amplitude of the mode (μ,ν) and τ is the MM-soliton pulse-width. In this case, the acquired angular frequency corresponding to each mode in this composite soliton is now denoted as \(\omega _{s_{\mu \nu }}\). In addition, the amplitudes \({\Phi}_0^{\mu \nu }\) can be determined from the following set of algebraic equations

$$\frac{{\left| {\beta _{\mu \nu }^{\left( 2 \right)}} \right|}}{{2\tau ^2}} = \mathop {\sum}\limits_{ {{\mu {\prime} \ne \mu } \atop {\nu {\prime} \ne \nu }} } {\left[ {R_{\mu \nu ,\mu \nu }\left| {{\Phi}_0^{\mu \nu }} \right|^2 + 2R_{\mu {\prime}\nu {\prime},\mu \nu }\left| {{\Phi}_0^{\mu {\prime}\nu {\prime}}} \right|^2} \right]} .$$
(4)

In obtaining the solution presented in Eq. (3), the hyperbolic secant ansatz was used, given that Eq. (2) is not integrable and hence cannot be in general addressed via inverse scattering transform methods. As opposed to SMFs, Eq. (4) clearly suggests that, for a given input energy, this MM-soliton solution is by no means unique. In other words, infinitely many modal distributions can be obtained, all satisfying Eq. (4). To some extent, all modes experience on average the same effective trapping nonlinear potential in a way akin to electron eigenfunctions in heavy atomic elements, when viewed within the framework of the Hartree-Fock method39. Once all the modes composing the MM-soliton are locked together at a common speed V, one can then determine the frequency shift Ωpq the mode (p,q) should undergo from the linear dispersion relation \(V^{ - 1} = \beta _{pq}^\prime ( {\omega _{s_{pq}}} ) = \beta _{pq}^\prime ( {\omega _0} ) + \beta _{pq}^{(2)}( {\omega _0} ){\Omega}_{pq}\), where \(\omega _{s_{pq}} = \omega _0 + {\Omega}_{pq}\) is the resulting carrier frequency of the participating mode (p,q), and ω0 denotes the carrier angular frequency of the input pulse. It is important to note that, while Eqs. (34) do represent an exact solution to this problem, in reality the temporal profiles of all modes involved in a MM-soliton could be considerably more complex. Yet, even in this case, once all modal components are locked together at a common group velocity V (because of mutual self-trapping), the corresponding frequency shifts Ωpq still remain the same as per our discussion above.

Starting from Eqs. (24), one can readily obtain the dispersion relationship at any angular frequency ω for each constituent spatial mode (p,q) involved in a MM-soliton:

$$\beta _{s_{pq}}\left( \omega \right) = \beta _{s_{pq}}\left( {\omega _{s_{pq}}} \right) + \frac{1}{V}\left( {\omega - \omega _{s_{pq}}} \right) + \frac{{\left| {\beta _{pq}^{\left( 2 \right)}} \right|}}{{2\tau ^2}},$$
(5)

where again \(\omega _{s_{pq}}\) represents the carrier frequency of mode (p,q) in the MM-soliton. The last term in Eq. (5) stands for the nonlinear shift in the propagation constant that each mode experiences during propagation. We note that in weakly guiding structures (like the ones considered here), the second-order dispersion does not appreciably change from mode to mode since it is dominated by material dispersion.

In general, dispersive waves result when a soliton sheds energy into a narrowband of frequencies through the interplay between higher-order dispersion effects and nonlinearity. In SMFs, the resonant frequencies of the emerging DWs are dictated by a phase-matching condition with respect to the propagation constant of the soliton itself. However, as we will see, unlike SMFs, the modal group delays in MMFs now play a deciding role in establishing the available phase-matching paths needed to generate DWs. Based on these latter arguments, the phase-matching condition between a soliton mode (p,q) and a DW in another mode (m,n) can only be met at a resonance frequency ωr when \(\beta _{s_{pq}}(\omega _{r_{mn}}) = \beta _{mn}( {\omega _{r_{mn}}})\), where βnm stands for the linear dispersion relationship associated with this mode. In practice, the soliton fission process further encourages the onset of DWs. For parabolic-index multimode fibers, the dispersion relation of the spatial mode (m,n) is given by \(\beta _{mn}\left( \omega \right) = k_0n\left( \omega \right)\left[ {1 - \left( {2\sqrt {2{\Delta}} (m + n + 1)/k_0an(\omega )} \right)} \right]^{1/2}\)where k0 = ω/c denotes the free space wavenumber at an angular frequency ω and n (ω) represents the refractive index of silica as function of frequency – as obtained from a corresponding Sellmeier equation. Here, a and Δ stand, respectively, for the core radius and relative index difference of the fiber. It is important to note, that, this latter dispersion relation was obtained within the Gauss-Hermite base of modal wavefunctions LPmn allowed in the parabolic-index MMF. This same system can also be considered in an equivalent Gauss-Laguerre modal base LPlp, provided that \(\left( {m + n} \right) \to \left( {2p + l - 2} \right)\). The group velocity of each mode in a graded-index potential can then be obtained from \(v_{g_{m,n}}^{ - 1} = d\beta _{mn}/d\omega\),

$$v_{g_{m,n}}^{ - 1} = \, \left[ {\frac{1}{c}\left( {n\left( \omega \right) + \omega \frac{{dn}}{{d\omega }}} \right)} \right]\left[ {1 - \frac{{\left( {m + n + 1} \right)\sqrt {2{\Delta}} }}{{\frac{\omega }{c}an\left( \omega \right)}}} \right]\\ \, \times \left[ {1 - \frac{{2\left( {m + n + 1} \right)\sqrt {2{\Delta}} }}{{\frac{\omega }{c}an\left( \omega \right)}}} \right]^{ - 1/2}.$$
(6)

Clearly, if the velocity V of a MM-soliton is known, then the frequency shift that is nonlinearly acquired by each of its constituent spatial modes (needed to compensate for the intermodal group-velocity mismatch), can be evaluated. On the other hand, the linear dispersion relation corresponding to a DW in mode (m,n) in such a parabolic-index MMF, can be approximated to first-order as \(\beta _{mn} = k_0n\left( \omega \right) - \left( {m + n + 1} \right)\sqrt {2{\Delta}/a^2}\). From here, the phase-matching condition between a DW in mode (m,n) and a component residing in mode (p,q) of this MM-soliton (traveling at velocity V) can be directly determined from (see Supplementary Note 1)

$$\omega _{r_{mn}}\left( {\frac{{n\left( {\omega _{r_{mn}}} \right)}}{c} - \frac{1}{V}} \right) = \left[ {\left( {m + n} \right) - \left( {p + q} \right)} \right] \, \sqrt {\frac{{2{\Delta}}}{{a^2}}} \\ + \frac{{\omega _{s_{pq}}}}{c}n\left( {\omega _{s_{pq}}} \right) - \frac{1}{V}\omega _{s_{pq}} + \frac{{\beta _{pq}^{\left( 2 \right)}}}{{2\tau ^2}},$$
(7)

where again \(\omega _{r_{mn}}\) is the angular frequency of the DW. Equation (7) plays a central role in our subsequent discussion.

Simulations and experiments

In order to verify our theoretical analysis, we performed numerical simulations using gUPPE techniques that globally take into account waveguiding, and nonlinear and higher-order dispersion effects in fused silica40. The beam spot-size assumed at the input was taken to be 15 µm so as to excite a multitude of LP0p Gauss-Laguerre modes. As shown in Fig. 1, the input pulse, after being injected into this MMF, undergoes a continuous temporal contraction along with a spectral expansion at a propagation distance of ~ 8 cm. At this point, MM-soliton fission (which is marked by a sudden expansion of the spectrum) takes place (Fig. 1a), and subsequently the propagating pulse disintegrates into a sequence of secondary multimode solitons (Fig. 1b). Right after fission, there is an emergence of broadband DW-comb lines, occupying the spectral region between 500 and 1000 nm. A closer look at the spectral evolution, reveals that within a short distance after the MM-soliton fission, additional DW lines start to appear in clusters, whose discrete frequencies can be again determined by the secondary MM-solitons inducing them, as expected from Eq. 7. Figure 1c depicts the evolution of the first emitted MM (secondary) soliton during the first few centimeters after its emission. This figure clearly shows that this first emitted soliton (having the highest peak intensity) undergoes significant Raman-induced red-shifting just few centimeters after the MM-soliton fission event. During this period, this decelerating MM-soliton experiences multiple temporal compressions, each of which is accompanied by a spectral expansion around the soliton spectrum and an emission of low intensity DW lines, close to the visible region. Our model (Eq. 7) now allows us to explain all the primary DW line-combs appearing in Fig. 1, right after the soliton fission point. More specifically, once the MM-soliton velocity is known (and hence the frequencies of the participating modes \(\omega _{s_{pq}}\)), one can then directly predict the spectral locations of the DW lines from Eq. 7. From our simulations, the MM-soliton’s velocity right after the fission event is estimated to be 2.049 × 108 ms−1 (ng = 1.4631) which is approximately to that expected from a soliton that is centered close to the pump wavelength. This is because the Raman-induced frequency red-shifting is still negligible at this stage. Figure 2 depicts the spectrum immediately after the soliton fission along with the positions of all DW lines (dotted lines), as predicted from our theoretical model (see Supplementary Note 2). As this figure suggests, there is an excellent agreement between the DW frequencies predicted by Eq. 7 and our numerical results. Note that, because the fiber has been excited on center, only axially symmetric modes (LP0p modes) carry energy and hence only modes with even (m + n) actively participate in reshaping the spectrum. This same analysis can now be repeated for subsequent temporal compressions experienced by the first MM-soliton (secondary soliton) after it has undergone a Raman red-shifting. At these later stages, the composite soliton velocity and wavelength have been altered, which, in turn, results in a new set of phase-matching conditions and hence different sets of DW spectral lines. The group velocity of this same soliton, at approximately 1.5 cm right after the MM-soliton fission event, is now estimated from simulations to be 2.0465 × 108 ms−1 (ng = 1.4649), which corresponds to a slower soliton at around 161.2 THz (1.8598 µm), because of substantial Raman redshifting. Our theoretical predictions for the DW lines associated with this redshifted wavelength are again in good agreement with the results obtained from simulations (see Supplementary Note 2). Note that our analysis has so far solely dealt with the first emitted MM-soliton. However, DW lines produced by the subsequent solitons can also be obtained following the same analytical procedure. In most cases however, the first emitted soliton carries most of the total energy and hence is responsible for the primary features of the DW content. In other words, while knowledge of MM-soliton velocities can facilitate our understanding of this process, the primary DW lines are generated right at the MM-soliton fission point, i.e., at the pump wavelength. Therefore, the group velocities of the various modes comprising the MM-soliton can be directly extracted from Eq. (6).

Fig. 1: Temporal and spectral evolution of a multimode optical pulse in a graded-index multimode fiber.
figure 1

a Spectral evolution of a 400-fs pulse (150 nJ) at 1550 nm after propagating in a 20 cm long graded-index multimode fiber having a core diameter of 62.5 µm, NA = 0.275 as obtained from a gUPPE simulation. b Temporal evolution of this same pulse, under the same conditions. c A closer look at multiple fissions and of the red-shifting of the first emerging soliton right after the primary soliton fission event. The horizontal axes of b and c are expressed in a moving coordinate frame T = t – z/v42.

Fig. 2: A comparison of theoretically predicted and numerically simulated dispersive wave spectrum in a graded-index multimode fiber.
figure 2

Dispersive wave (DW) spectrum as obtained from a gUPPE simulation right after the first fission event. The parameters used are identical to those of Fig. 1. Dotted lines indicate the discrete positions of the DW lines, as predicted by Eq. 7.

To corroborate our theory, we performed a series of experiments with both parabolic and step-index fibers. The parabolic MMF used had a core radius of 35 µm and Δ~0.021 (NA ≈ 0.2). The fiber was illuminated on-axis with a 100 fs pulse (500 kW) at 1550 nm emitted from an OPO (1 kHz) pumped by a Ti:Sapphire laser. Figure 3 shows the spectrum recorded at the end of a 10 m long fiber along with the theoretically calculated DW line frequencies (dotted lines) from Eq. 7. For this analysis, the soliton velocity was computed from the Sellmeier series, by assuming that most of the MM-soliton spectral content resides around the pump wavelength. Because of the radially symmetric fiber excitation used, preservation of optical angular momentum demands that during conversion, the spectral features are dominated by DWs of even order, i.e., \(\left( {\left( {m + n} \right) - \left( {p + q} \right) = 2K,K = 0, \pm 1, \pm 2, \ldots } \right)\). Figure 3 reveals, that, from the fiber parameters, our theory can accurately predict the DW lines emerging right at the MM-soliton fission point (see Supplementary Note 3). Interestingly, in a same vein, the location of subsequent DW line combs emitted by secondary MM-solitons can spectroscopically provide information as to their respective velocities.

Fig. 3: A comparison of theoretically predicted and experimentally observed dispersive wave spectrum in a graded-index multimode fiber.
figure 3

Experimentally measured spectra at the output of a 10 m long parabolic-index multimode fiber of core radius 35 µm. Theoretically predicted dispersive wave wavelengths (Eq. 7) are shown in dotted lines.

According to Eq. (7), to populate the far reaching DW lines, the beam at the soliton fission point should involve a significant portion of the total energy in higher-order modes. Supplementary Fig. 2, obtained from the numerical gUPPE simulations, shows that at this point the beam is significantly compressed in both space and time, indicating that higher-order modes are indeed populated at soliton fission, irrespective of the modal composition at the input (see Supplementary Note 4). Given that the energy initially resides in the lowest order group of modes, our model predicts that DW lines lying in higher-frequencies must reside in higher-order modes in order to satisfy the phase-matching condition of Eq. (7). This is especially necessary, in exciting high frequency DW Cherenkov waves and it is in agreement with the experimental results reported in ref. 18. Similarly, one can also alter the spectral content of the DW lines by judiciously modifying the input modal composition20. In other words, by initially populating higher-order LP0p, the main portion of the DW energy can be transferred to higher frequency regions20.

Another interesting aspect anticipated by both our theoretical model and simulations is the appearance of a spectral DW line in the long-wavelength region (at 2.5 µm in Fig. 1) above the pump wavelength. Supplementary Fig. 3 depicts the phase-matching conditions for the fundamental mode, suggesting that they can be met in both the anomalous as well as in the normal dispersive regions (see Supplementary Note 5).

Both our experimental observations and numerical results indicate that the wideband generation of dispersive waves in graded-index MMFs always occurs in a discrete fashion—a direct consequence of the equidistant distribution of propagation eigenvalues. Equation 7 dictates that the phase-matching condition can be met only at discrete frequencies since the term (m+n)−(p+q) is always an integer. Since this feature does not apply for other index distributions, one may then expect a drift for each DW frequency and an absence of any clustering of lines. To verify this hypothesis, we simulated a step-index fiber (105 µm in core diameter and NA = 0.275). A similar fiber was also used in a different set of experiments. In Fig. 4, the visible portion of the experimentally measured spectrum for this step-index MMF (100 fs, 700 kW, 1550 nm, 10 m) is contrasted to that previously obtained in the parabolic MMF. Our experiments clearly show that DW generation plays out distinctly even in step-index fibers, an aspect that previous models could not account for. The numerically simulated spectrum, resulting shortly after the soliton fission event is also depicted in Supplementary Fig. 4 in Supplementary Note 6. These figures clearly indicate that indeed the DW grouping vanishes in step-index MMFs. As a result, the DW spectrum forms a continuum extending from 450 to 1000 nm. These results suggest that the distribution of the generated DWs in MMFs can also be tailored through appropriately designed refractive index profiles.

Fig. 4: A comparison between experimentally observed dispersive spectrum in a graded-index and a step-index multimode fiber.
figure 4

A comparison between dispersive waves (DWs) generated in a step-index multimode fiber (MMF) (blue curve) and a parabolic-index MMF (green curve). The step-index MMF had a length of 10 m, a core diameter of 105 µm and a numerical aperture NA = 0.22. In both cases the extent of the DW spectrum is approximately the same. However, the DWs generated from a step-index MMF do not have the comb-like features, a signature of parabolic MMFs.

So far, in our experiments, we have only considered moderate input power levels. As we apply larger input powers, the nonlinear phase term (Eq. 4) starts to play a more important role in the phase-matching condition (last term on the right-hand side of Eq. 7). To experimentally observe the impact of this nonlinear phase-contribution, in both graded-index and step-index MMFs, we gradually increased the input power injected and recorded the output spectrum. As it is shown in Fig. 5, because of this nonlinear phase term, the DW lines gradually drift toward higher frequencies as analytically predicted by Eq. 7. However, since lower order modes have larger nonlinear coefficients, their spectral shift happens to be more significant.

Fig. 5: Experimental results indicating the blue-shifting of dispersive waves resulting from the nonlinear phase term.
figure 5

a Experimentally measured blue-drift of dispersive waves (DWs) in a parabolic multimode fiber (MMF) as a result of increasing the input power levels. The blue drifting of DWs arises from the nonlinear phase term (last term in Eq. 7). b Experimentally measured blue-drifting of DWs for a step-index MMF when the input power increases. In both cases, a similar trend is observed. We note that the maximum power level used in these experiments 1.7 MW is well below, the threshold power required at 1.55 µm for catastrophic self-focusing (10 MW) in silica fibers.

In conclusion, we have provided a comprehensive theory capable of explaining the distinct DW (Cherenkov) lines emerging after multimode soliton fission events. Our results indicate that this broadband DW emission results from the nonlinear self-trapping of the constituent modes involved in the generated multimode solitons. In this respect, DWs can be generated in both the normal and anomalous dispersion regions of a MMF. Experiments carried out in both parabolic and step-index multimode optical fibers are in good agreement with our theoretical model. Our results could pave the way towards understanding the complex dynamics of highly multimode and multidimensional nonlinear systems.

Methods

Simulations

To verify our theory, we conducted simulations using generalized unidirectional pulse propagation equations (gUPPE). For the case of graded-index fiber, the core radius of the fiber was assumed to be 31.25 µm, the numerical aperture was assumed to be NA = 0.275 and the silica Kerr coefficient was taken to be 2.9 × 10−20 m2W−1. For the case of step-index fiber, the core diameter of the fiber was assumed to be 105 µm, the numerical aperture was assumed to be NA = 0.22 and the silica Kerr coefficient was taken to be 2.9 × 10−20 m2W−1. In addition, all nonlinear and linear processes such as modal walk-offs, self-focusing, self-phase modulation, cross-phase modulation, FWM, Raman and shock effects were accounted for in the simulations. The input pulse in the anomalous regime was centered around 1550 nm (0.193 PHz). In all our simulation, the pulse-widths were assumed to be 400 fs and the longitudinal step size was set at Δz ≈ 1 µm. For undertaking such computationally demanding simulations, we used parallel computing, provided by the XSEDE supercomputer facility.

Experimental design

To corroborate our theory, we performed a series of experiments with both parabolic and step-index fibers. The parabolic MMF used was 10 m long, had a core radius of 35 µm, and Δ ~ 0.021 (NA ≈ 0.2). The fiber was illuminated on-axis with a 100 fs pulse at 1550 nm emitted from an OPO (1 kHz) pumped by a Ti:Sapphire laser. The step-index fiber used in our experiments was in-house fabricated germanium-doped silica multimode optical fiber. The fiber has a core diameter of 105 µm and a length of 10 m. The 100 fs pulses were coupled to the from a Ti:Sapphire laser source (1.55 μm). A three-axis translation stage and a 5 cm focal length lens were used to couple laser beams into the fibers. The ensued spectrum was measured with visible and IR spectrometers. In acquiring our data, both polarizations are accounted for.