General theory and observation of Cherenkov radiation induced by multimode solitons

Advancements in computational capabilities along with the possibility of accessing high power levels have stimulated a reconsideration of multimode fibers. Multimode fibers are nowadays intensely pursued in terms of addressing longstanding issues related to information bandwidth and implementing new classes of high-power laser sources. In addition, the multifaceted nature of this platform, arising from the complexity associated with hundreds and thousands of interacting modes, has provided a fertile ground for observing novel physical effects. However, this same complexity has introduced a formidable challenge in understanding these newly emerging physical phenomena. Here, we provide a comprehensive theory capable of explaining the distinct Cherenkov radiation lines produced during multimode soliton fission events taking place in nonlinear multimode optical fibers. Our analysis reveals that this broadband dispersive wave emission is a direct byproduct of the nonlinear merging of the constituent modes comprising the resulting multimode soliton entities, and is possible in both the normal and anomalous dispersive regions. These theoretical predictions are experimentally and numerically corroborated in both parabolic and step-index multimode silica waveguides. Effects arising from different soliton modal compositions can also be accounted for, using this model. At a more fundamental level, our results are expected to further facilitate our understanding of the underlying physics associated with these complex “many-body” nonlinear processes. The multifaceted nature of nonlinear multimode optical fibers has provided a fertile ground for observing novel physical effects otherwise impossible in single-mode setting. The authors present a theoretical and experimental study, providing a comprehensive theory capable of explaining the distinct nature of wideband Cherenkov radiation produced during multimode soliton fission events.

T he past few years have witnessed a resurgence of interest in multimode waveguide structures, predominantly driven by the ever-increasing demand for higher information capacities [1][2][3][4] . This renaissance, in turn, incited a flurry of activities in the general area of nonlinear multimode fiber optics  . 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 settings 12 . These include for example, the formation of multimode solitons (MMsolitons) 11 , beam self-cleaning effects 16,17,19 , geometric spatiotemporal instabilities 14 , and efficient supercontinuum and second-harmonic generation 17,20,21,25 , to mention a few. The prospect for spatiotemporal laser mode-locking has also been successfully demonstrated in a recent study 26 .
One of the most prominent processes taking place in nonlinear optical waveguides is that associated with dispersive wave (DW) emission or Cherenkov radiation [27][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 microstructures [33][34][35] .
As one could expect, this situation is considerably more involved in nonlinear multimode environments. In fact, in a recent experimental study 12 , 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 solitons 36 . 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 selfimaging 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 r; z; ω ð Þ¼n 0 r; ω ð Þþn 2Ê r; z; t ð Þ 2 , where n 0 (r, ω) represents the waveguide index profile and n 2 the nonlinear Kerr coefficient of silica glass. Following the above arguments, we express the total electric fieldÊ r; z; t ð Þ in this system as a superposition of all the spatial modes involved, i.e., 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 Porto 37,38 , one then finds that, in the anomalous dispersion region, the envelopes evolve according to represent the second-order dispersion coefficient and modal-group velocity. The self-phase and cross-phase modulation coefficients are given by R μ0ν0; 38 . In developing Eq. (2), we intentionally omitted any four-wave mixing (FWM) terms. This FWM suppression can be justified 37,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, where Φ μν 0 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 ω s μν . In addition, the amplitudes Φ μν 0 can be ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00640-1 determined from the following set of algebraic equations 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 method 39 . 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 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. (3)(4) 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. (2)(3)(4), one can readily obtain the dispersion relationship at any angular frequency ω for each constituent spatial mode (p,q) involved in a MM-soliton: where again ω 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 β s pq ðω r mn Þ ¼ β mn ðω 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 k 0 = ω/c denotes the free space wavenumber at an angular frequency ω and n (ω) represents the refractive index of silica as function of frequencyas 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 LP mn allowed in the parabolic-index MMF. This same system can also be considered in an equivalent Gauss-Laguerre modal base LP lp , provided that m þ n ð Þ! 2p þ l À 2 À Á . The group velocity of each mode in a graded-index potential can then be obtained from v À1 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 . From here, the phasematching 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) where again ω 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 silica 40 . The beam spot-size assumed at the input was taken to be 15 µm so as to excite a multitude of LP 0p 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 MMsoliton 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 redshifting 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 ω 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 × 10 8 ms −1 (n g = 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 (LP 0p 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 × 10 8 ms −1 (n g = 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 MMsoliton 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).
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., m þ n ð ÞÀ p þ q À Á ¼ 2K; K ¼ 0; ± 1; ± 2; À Á . 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.
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 composition 20 . In other words, by initially populating higher-order LP 0p , the main portion of the DW energy can be transferred to higher frequency regions 20 .
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.
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 stepindex 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.
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 stepindex 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 m 2 W −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 m 2 W −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.

Data availability
All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.