Spatiotemporal dynamics of multimode optical solitons

We study multimode optical solitons with up to roughly 10 spatial modes. This work provides the first evidence for solitons involving more than a few modes, and for spatiotemporal multimode soliton fission and Raman shifting.


Background
Single-mode optical fibers (SMFs) provide a convenient platform for the study of the solitons of the nonlinear Schrödinger equation (NLSE), which are now among the most thoroughly studied phenomena in nonlinear dynamics [1,2]. The formation and interaction of these objects provides a conceptual framework for understanding many nonlinear optical phenomena, from supercontinuum generation in SMFs [3,4] to mode-locked lasers [5,6] and optical rogue waves [7]. In multimode systems like multimode fibers (MMFs), more complex multimode solitons (MM solitons) are theoretically possible [8][9][10][11][12][13][14][15][16]. For the purposes of this article, we will use the following definition: MM solitons are non-dispersive pulses whose energy is distributed over multiple spatial modes. They hence counter both intra-and intermodal dispersion (i.e., group velocity dispersion and modal velocity mismatch) and should come in a variety of different spatiotemporal shapes, depending on the modes involved. Like (1+1)-D solitons for SMF, MM solitons may provide a simple framework for the interpretation of complex nonlinear dynamics in MMFs. This is a key reason for interest in MM solitons. They should exhibit multimode spatiotemporal dynamics unlike single-mode solitons, and represent an important potential platform for applications in multimode fiber lasers and spatial division multiplexing.
In spite of considerable theoretical work [8][9][10][11][12][13][14][15][16][17], to date there is only one experimental study [18] of MM solitons in MMFs. Other studies on nonlinear optics in MMFs have focused primarily on few-mode intermodal frequency mixing effects, as modal dispersion allows phase-matching of efficient, broadband four-wave mixing [19][20][21] and dispersive wave radiation [22][23][24][25]. Studies in the normal dispersion regime have shown strong intermodal Raman effects [26][27][28]. Intermodal nonlinear effects are also important in gas-filled photonic crystal fibers, where Raman-free tunable nonlinearity, tunable dispersion and plasma effects make for a powerful experimental platform [24,25]. The recent work by Renninger et al. [18] studied solitons in a graded-index (GRIN) MM fiber, whose relatively small intermodal dispersion facilitates intermodal nonlinear effects. This study provided the first experimental evidence towards confirming the numerous theoretical predictions [8,9,[11][12][13]17] of and about multimode solitons. However, as the soliton was comprised primarily (>90%) of the fundamental mode and contained a relatively small total energy (0.5 nJ), nonlinear effects -including the synchronization of the pulses from two higher spatial modes -were due almost exclusively to the fundamental mode. Consequently, the soliton behaved very similarly to a single-mode soliton. Many properties and behaviors of the full range of MM solitons hence remain unknown or unconfirmed by experiment. First, it will be important to confirm that solitons with higher energy and with higher energy in higher-order modes (HOMs) exist. Second, we should investigate how MM solitons form from a given launched pulse, and understand the analogue of soliton fission. Similar to the fission of polarization vector solitons [29], we should expect the soliton fission process in MMFs to result in MM solitons and dispersive radiation with various spa-tiotemporal distributions. Finally, the response of MM solitons to perturbative effects such as Raman scattering should be studied. These aspects will help to build up a knowledge base that may eventually allow for intuitive interpretations of nonlinear multimode dynamics such as supercontinuum in MM fibers. In this work, we study numerically and experimentally these issues in a multimode GRIN fiber and provide a first look at the complex spatiotemporal dynamics of MM solitons. In doing so, we also provide promising indications for the application value of multimode optical solitons.

Methods
We use two complementary approaches to numerically model propagation of spatiotemporal optical pulses in the multimode GRIN fiber. The first and most direct is a generalized (3+1)-D nonlinear Schrödinger equation (GNLSE): = n e f f ω/c, n e f f is the average index,D is the Taylor series approximation to the chromatic dispersion operator = h R is the Raman response of the fiber medium, f R ≈ 0.18 is the Raman contribution to the Kerr effect , γ = n 2 k e (ω o ) and n 2 is the nonlinear index. In the case of a quadratic index profile (and neglecting Raman, shock and dispersion orders greater than β 2 ), Eq. (1) is a Gross-Pitaevskii equation, which is used to model Bose-Einstein condensates [30]. Alternatively, one may employ a set of coupled (1+1)D GNLSEs, for N modes of consideration, each with temporal envelope A p . This is known as the generalized multimode NLSE (GMM-NLSE) [31, 32]: , where δ β 1 ) are the differences relative to the fundamental mode of the 1st (2nd) coefficient in the Taylor series expansion of the wavenumber for the mode p about ω o ,D p iŝ D as defined above for mode p, and S R plmn and S K plmn are the nonlinear coupling coefficients for the Raman and Kerr effect respectively [32]. For all simulations using Eq. (2) discussed in the present work, we assume all modes share a single linear polarization and neglect the shock term. For all simulations we use n 2 = 2.3 × 10 −20 m 2 /W , and in Eq. (1) use β 2 = -25 f s 2 /mm and β 3 = 140 f s 3 /mm. For Eq. (2), the spatial modes and dispersion parameters are calculated at 1550nm following the same method as in previous work [18], modelling the chromatic dispersion of the core index by Ge-doped fused silica [33]. The GMM-NLSE has, for example, been used to model supercontinuum, self-focusing and FWM in MMFs [33-35].

Example simulation results
Before introducing our experimental results, we consider numerical simulations of pulse propagation using Eq. (2). These simulations show qualitatively the different energy regimes and characteristic nonlinear dynamics in each ( Fig. 1) and will serve as a conceptual guide throughout the rest of the article. We consider for simplicity an arbitrary particular initial condition   where the first five modes of the GRIN fiber are excited with identical pulses with equal energy. As will be discussed below, details of the observed behavior depend on the initial launch conditions (distribution among modes). For example: the modal distribution of the soliton(s) formed and the energy level at which fission first occurs both depend on the initial mode distribution. However, for nearly all conditions we have considered, the main features of the observed dynamics are qualitatively similar. We consider only the modes without angular momentum (i.e., m = 0 in the notation used by Mafi [33]). These modes are the ones excited by a centered Gaussian beam, and over short distances are relatively unaffected by linear mode coupling from, e.g. bending [33]. Simulations conducted with a mixture of m = 0 and m = 0 modes produced qualitatively similar results as those with only m = 0 modes. We assume a GRIN fiber with 0.275 numerical aperture (NA) and a 62.5 μm core diameter. We launch 1550 nm, 500-fs Gaussian pulses and monitor propagation through 25 m of fiber. At very low pulse energy, pulses in each mode disperse individually due to intramodal (group velocity) dispersion and separate from one another due to intermodal dispersion ( Fig. 1, top row). As pulse energy increases, a single MM soliton may form, consisting of synchronized pulses from several modes. The pulse in each mode synchronizes due to wavelength shifts For simulated spectra, the full-field spatially-averaged spectra are shown, which accurately represents the spectral energy distribution in the field. The energies are 3.2, 5, 10, 14, 24 and 34 nJ (experiment) and 0 (i.e., linear propagation), 1, 10, 14, 19 and 22 nJ (simulation). The purpose of the comparison is to show the qualitative similarity; the two plots correspond to different initial modal energy distributions, the experimental launched pulse has a non-Gaussian spectrum, and there is a strong wavelength-dependent loss which is not considered in the simulation. Nonetheless, the quantitative agreement is acceptable until loss becomes important (when the Raman soliton exceeds 1800 nm). The energy in the most red-shifted peak is calculated to vary from 0 to approximately 40%, in agreement with a measurement made using a Ge window. caused by intermodal nonlinear effects. In the GRIN fiber considered here, lower-order modes shift to redder wavelengths in order to speed up and reach a common group velocity with pulses in HOMs. Note that the cause of this wavelength shift can be intermodal Raman scattering, or purely Kerr effects such as intermode cross-phase modulation or four-wave mixing (as in Ref. 18, where Raman was neglected). In Fig. 1 (middle row, energy regime 1: roughly < 5 nJ), most of the soliton's energy is within the first 4 modes. At greater energy, the pulse undergoes fission into multiple MM solitons and MM dispersive radiation ( Fig. 1 bottom row, energy regime 2: roughly > 5 nJ). The resulting MM solitons show a remarkable spatiotemporal diversity: many different combinations of modes can comprise a MM soliton. The most intense soliton typically contains more of the fundamental mode than any higher-order mode (HOM), but also usually contains non-negligible energy from HOMs -it remains a MM soliton. In the fission process and subsequently, stimulated Raman scattering leads to amplification of the lower-order components of these MM solitons. Their spectrum shifts red due to a process resembling the well-known soliton self-frequency shift (SSFS) [36, 37], but pulses in all the modes shift their spectrum in a staggered manner, preserving their temporal synchronization. The Raman scattering bias towards the fundamental mode is due to two cooperating factors. First, because the low-order modes in the MM solitons are redshifted to form the soliton, they are spectrally positioned to receive Raman gain from the bluer HOMs. Secondly, of all the modes in a GRIN fiber, the fundamental mode has the most overlap with the other modes [38]. This typically causes it to experience the highest Raman gain in a multimode amplifier, an effect commonly called 'Raman beam cleanup' [38,39]. Together, these factors account for the strong intermodal energy transfer, via Raman scattering, into the fundamental mode. The combined effect produces an ultrashort, spatiotemporally-localized pulse red-shifted from the pump field and was evidently observed in early experiments in the >1300 nm emission produced by excitation with 90 μJ, 150 ps, 1064 nm pulses. However, the authors were unable to fully explain their observations [27].

Experimental methods
Experimentally, we launch 500-fs pulses at 1550 nm from an amplified erbium fiber laser into GRIN fibers with 0.275 NA and 62.5 μm core diameter. A lens, half-wave plate, polarizing beamsplitter, and three-axis translation stage allow control of the initial energy and modal excitation. The important distinctions from previous work are that a) the fiber has a graded-index profile, b) we deliberately attempt to excite combinations of, rather than individual, modes, and c) we interpret our results in the context of simulations which show clear multimode soliton behavior. Mode coupling due to mechanical waveguide perturbations is minimized by primarily centered excitation (exciting m = 0 modes), by use of a relatively short, jacketed, connectorized fiber, and by minimizing bending [33]. As the spectrum of the field varies with space, a given measurement of the spectrum depends on the alignment and may not reflect accurately the total spectral energy distribution. Therefore, we record an overall averaged spectrum by using a defocused lens to couple the light into single mode fiber before measuring the spectrum. By using a germanium window, which acts as a long-pass filter cutting on at roughly 1800 nm, we can determine the amount of energy that is shifted beyond 1800 nm. These two spectral measurements (the spectrum analyzer, and the power meter after a long-pass filter) agree with one another and with simulations, indicating that the Raman soliton at wavelengths greater than 1800 nm contains 10-50% (depending on launch conditions) of the total launched energy. To measure the pulse in the time domain, use of a two-photon autocorrelator is critical as it avoids the problem of phase-matching a complex, broadband spatiotemporal pulsed beam. Figure 2 (left panel) shows the spectrum of the field generated by solving Eq. (2) with the initial condition mentioned above, for increasing energy. Figure 2 (right panel) shows the spectra measured for increasing energy in an experiment with a 25-m GRIN fiber. As the launch conditions and initial pulse spectra are by no means identical and because the simulation neglects infrared-wavelength loss, the comparison is meant only to show the qualitative similarity. Neglecting the energy scale (in Fig. 2, the experimental 2100-nm Raman soliton has 14 nJ energy and 90 fs duration), both plots are also qualitatively similar to the result which would be expected in a single-mode fiber [36,40]. It is natural to wonder if the experimental results can be explained by propagation of the fundamental mode only. Pulses in HOMs excited initially could simply disperse, rather than becoming involved in a MM soliton. In this case, the beam at the output of a fiber would still show multiple modes. However, the modes would not be nonlinearly coupled.
In the following sections, we consider an experimental test to confirm multimode soliton behavior and to disprove the null hypothesis just mentioned. The results of this test indicate that multimode solitons are observed, with similar behavior as expected from simulations. Further analysis of the experimental results provides some deeper insights into the behavior of multimode solitons. Then, we turn to a higher energy regime and study the spatiotemporal effects of multimode soliton fission and Raman scattering. The present work is the first to seriously consider multimode pulse propagation and multimode solitons except under extremely limited and simplified conditions. Furthermore, direct spatiotemporal or modal measurements (e.g., [41][42][43][44][45][46][47][48][49]) are not practical given the relatively low energy and average power, the relatively large number of modes, and the relatively poor performance of InGaAs technology in comparison with Si used for shorter wavelengths. Therefore, in order to interpret the results correctly, we need to account for and verify several different assumptions. While this analysis is necessary to justify the conclusions, the conclusions themselves can be understood without consideration of the detailed methods and analysis. These conclusions have been summarized in the first paragraph of most sections below.

Outline of a test to experimentally observe multimode solitons
In this section we outline a test to prove that multimode solitons are formed. This test is conducted in energy regime 1 (Fig. 1 middle row), where a single MM soliton forms from the launched pulse. The core of our test is that the required nonlinear phase shift needed to balance dispersive effects and form a MM soliton of a given duration should be at least equal to the energy required for a single-mode soliton of equivalent spacetime localization. This test confirms MM solitons without need for simulations. However, for an added layer of evidence, we also perform and find excellent agreement with a similar numerical experiment using Eq. (1). As in previous works, we consider only one transverse dimension [11][12][13]18] and assume a radial symmetry of the results in order to compare with experiment. The main conclusion of this section is that MM solitons are observed in our experiments, as in simulations, for a wide range of sizes and mode distributions, including up to 8-13 modes and with average spatial dimensions up to ≈ 3 times the fundamental mode size.
In the following analysis two implicit approximations are made. At the end, we will confirm that they are both justified. The first is that the experimentally measured beam size is roughly the average beam size, which in general a MM soliton will oscillate around. For the fiber considered here, this oscillation has a period of about 400 μm [18]. The phase of this oscillation depends on many conditions of the experiment which fluctuate over the many pulses that contribute to an experimental measurement (e.g., the phase front and position of the beam as it enters the fiber, or as it exits the fiber), and may vary across the wavelength of the broadband pulse. The phase also depends on conditions which vary uncontrollably from experiment to experiment (e.g., the fiber length to within 100 μm or the precise launch position). Hence over the many pulses that contribute to most measurements, and over many experiments, spatial measurements should tend to reflect the average spatial dimensions. Second, we assume that, provided the pulse does not split apart (which occurs in a higher energy regime discussed in a later section) the total energy is relatively close to the energy contained within the soliton. Since the pulse will emit dispersive radiation as it reshapes into a soliton, this approximation is not completely accurate. Nonetheless, it will turn out to be accurate enough for our present purpose of proving the observation of MM solitons and investigating their qualitative behavior.
For a single-mode soliton, the relationship between the energy, E p , and duration, τ, is: where λ o is the center wavelength. R g is the spatial dimension (waist) of a Gaussian beam or mode, so that E p /(R g τ) is the pulse intensity. For beams in a multimode fiber, different combinations of modes can produce a total beam with R g that is on average greater than or equal to the fundamental mode size, R o . To form a MM soliton, nonlinearity must somehow balance intramodal dispersion (group velocity dispersion) as well as intermodal dispersion (modal velocity mismatch). The nonlinear phase shift φ ∝ n 2 I ∝ n 2 E p /(τR 2 g ) should therefore be greater than or equal to the nonlinear phase shift required for a single-mode soliton of the same spatiotemporal localization τR 2 g . For a fixed R g , a plot of measurements 1/τ as a function energy, E p , may be fitted by a straight line. The line's slope, M, characterizes the relationship between the nonlinear and linear dispersive phases at balance. If the observed pulses are MM solitons, then they will require relatively more nonlinear phase for a given τ. Hence the single-mode soliton slope, M 1 (R g ) (Eq. (3)), should be a lower bound for the slope M(R g ) measured for a MM soliton.
As a null hypothesis, consider the hypothetical scenario in which the only observed solitons are purely in the fundamental mode, and pulses in all other modes only broaden and move away (i.e., propagate linearly). That is, the field is a fundamental mode soliton plus a higher-order mode background (HOMB). This is an important case to consider, as measurements would imply R g > R o , and M > |β 2 |R 2 o λ o /n 2 (i.e. the slope for a single-mode soliton in the fundamental mode only), but no MM solitons would be present. For this situation it is straightforward (see Appendix) to derive the measured slope as modified by the effect of the HOMB, M FHB , as a function of R g . However, the precise form of M FHB (R g ) is unimportant. The important conclusion is simple: since the energy within the actual single-mode soliton is less than the total measured energy, the slope measured for a single-mode soliton plus HOMB, M FHB (R g ), is always less than or equal to the single-mode soliton slope M 1 (R g ). Therefore, as long as the measured slope M(R g ) exceeds M 1 (R g ), it automatically exceeds M FHB (R g ). Furthermore, by additional analysis (see Appendix) we can show that the presence of a HOMB will cause the slope measured for a MM soliton to decrease, and therefore should only cause a false negative rather than a false positive result. This proves that the presence of a HOMB does not invalidate the test.
In experiments and simulations, 500-fs, 1550-nm pulses are launched into GRIN fibers as Gaussian beams of various sizes. These excite MM solitons with various sizes R g . The nearfield beam profile is measured and fit with a Gaussian. For each R g , the energy is varied, and for each value of E p the pulse duration is measured. Note that, for the energies considered in this section (energy regime 1 in Fig. 1), changing E p did not cause R g to change in any significant or consistent way. The E p vs 1/τ data is fit by a line to retrieve the slope M(R g ).
For the simulations, in order to account for the possible spacetime coupling of the field, we compute the spatially-integrated interferometric autocorrelation which for an uncoupled field simplifies to the usual result. In all cases, a Gaussian fit is used to characterize the real spatial dimension of the field. As the fit is only used to determine the equivalent-size Gaussian beam, the beam profile needn't be precisely Gaussian. However for most cases the observed beam profiles are relatively Gaussian, reflecting the center launch condition. For the simulation, R g is the average over the last ≈ 50 cm of fiber.

Results of the test to observe multimode solitons
When this test is carried out, we find that the experimental results confirm observation of MM solitons, and furthermore agree quite well with the simulation results. Mathematically, the slope of E p versus 1/τ, M, should fulfil the condition M(R g ) ≥ M 1 (R g ). Figure 3 shows the resulting M vs R g values for both simulation and experiment, as well as the expected relationship for single-mode solitons, M 1 (R g ) (Eq. (3) to the spatiotemporal one used so far. In Fig. 3, the simulation and experimental points are overall similar in their trend and magnitudes. Since the simulation values of R g were obtained by averaging over the beam oscillations, the assumption that the measured beam size corresponds to the average beam dimension is supported. In both cases the measured beam size could be inflated by dispersive radiation in HOMs, so the value may actually slightly exceed the average size of the soliton. However, since this would only increase the value of M required to prove MM soliton behavior, the test results remain valid. As is apparent in Fig. 3 and from Eqs. (8) and (9), energy in a HOMB causes the measured value of M for either a single-mode soliton or MM soliton to increase quite slowly with R g compared to the observed behavior for both MM and single-mode solitons. In essence, the test M(R g ) ≥ M 1 (R g ) is sufficient to prove MM solitons, whether a HOMB of dispersive radiation is present or not. Therefore, neglecting the dispersive HOMB is justified.
The beam size of the MM solitons oscillates and so a certain modal distribution corresponds to a range of beam sizes. Figure 4 shows the number of modes contributing more than 0.5% of the total energy to a multimode beam with a Gaussian spatial profile as a function of R g . The points are extracted by computing the mode overlap integrals for the first 30 m = 0 modes of the GRIN fiber with a Gaussian beam with varying R g . There are two branches of points, which correspond to the maximum and minimum beam sizes. To extract the average beam size as a The y-axis shows the number of modes that contribute more than 0.5% of the total energy the beam with size R g . To provide a continuous measure of the mode content as a function of average beam size, we fit the points, and then averaged the resulting curves to produce the average curve. For Gaussian-shaped beams/solitons, this line represents the number of modes contributing as a function of the average spatial localization of the beam/soliton. function of mode content, we fit curves to each set, and then average these curves. The resulting average curve provides a measure of the number of modes required to represent a beam of a given average spatial localization. Since we are presently considering only m = 0 modes, a practical upper limit on the number of modes in a MM soliton is approximately 20, as beyond that the maximum spatial dimension would exceed the core dimensions. The maximally-multimode solitons in the present work contain between 8-13 modes. More importantly, the MM solitons exhibit characteristics that are quite distinct from solitons dominated by the fundamental mode [18].

Fission and Raman scattering: energy regime 2
Having now explored the first energy regime, where only a single MM soliton forms and where higher-order effects are not important, we turn to the higher energy regime. This regime, corresponding to the bottom row of Fig. 1, emphasizes the spatiotemporal aspects of MM solitons. For the range of conditions considered here, the typical sequence is as follows. The launched MM pulse first compresses, then undergoes fission into numerous MM solitons, each containing a different mode distribution, as well as MM dispersive radiation. The most intense soliton is redshifted by intrapulse Raman scattering, while preserving its multimode characteristic. The transition from the previous energy regime occurs at a total energy that evidently depends on the initial mode energy distribution, and it is apparent experimentally as the spectrum shows the formation of a clear red-shifting peak. Furthermore, additional pulses are usually observed and the duration of the main peak of the autocorrelation ceases to decrease with increasing energy, or decreases at a slower rate than observed for lower energy launch pulses. This behavior is qualitatively similar to single-mode fiber solitons. However, for MM solitons the processes are spatiotemporal and, owing to the diversity of possible MM solitons, more varied and rich. Simulations based on Eq. (1) with one transverse dimension illustrate this behavior in 2 m of GRIN fiber (Fig. 5). The left (middle) panel shows the expected interferometric autocorrelation for several pulse energies when the initial launched beam has a 10 μm (15 μm) waist. The autocorrelation of a spatiotemporal field (Eq. (4)) has some limited spatial information. Particularly when the Raman effect is prevalent, the relatively higher intensity (lower-order mode content) of the Raman soliton as compared with the remnant components causes the secondary peaks observed in the autocorrelation to decrease in amplitude. The degree to which secondary peaks are reduced relative to their actual energy content reflects how much less spatially localized they are than the Raman MM soliton. The right panel shows the spatiotemporal spectrum corresponding to the uppermost autocorrelation in the middle panel (28 nJ). After 4 cm, the spectrum shows that numerous modes are excited, evidenced by the non-Gaussian angular distribution (top). However, the spectral and spatial distributions are still separable. Following soliton formation and fission, the redshifted portion of the spectrum, which corresponds to the accelerated Raman soliton, exhibits relatively little structure as compared to the remnant MM solitons and non-solitonic radiation. Evidently, the field is spatiotemporally coupled and complex (Fig. 5 right bottom).
Our experimental measurements show similar spatiotemporal phenomena. We are spectrally limited by the 1750-nm absorption edge of the InGaAs camera used to measure the beam profile. Hence, in order to observe a strong spatial nonlinear effect, we need significant initial HOM content in order to ensure the primary elements of the dynamics occur within the spectral range 1550 nm to 1750 nm. Therefore the alignment into the fiber is adjusted until the lowenergy output beam profile has significant energy in HOMs (two-peaked and elliptical, Fig. 6). The experimental measurements of the spectrum, beam profile and interferometric autocorrelation, as well as their trends with increasing energy, all agree qualitatively very well with the simulations. Consequently, following the simulations we interpret the dynamics as follows. At 4 nJ, the interferometric autocorrelation (Fig. 6 left) suggests a MM soliton is formed (with secondary features in the AC corresponding to HOMB pulses, as in the middle row of Fig.  1). As the energy is increased to 7 nJ, a larger fraction of the launch energy is coupled into the MM soliton, and the pulse compresses. At higher energies, the pulse undergoes fission, and the most intense soliton produced is Raman shifted. Secondary peaks in the autocorrela-tion are quite weak, as in the simulation results. These weaker, localized pulses correspond to less-energetic MM solitons (or those containing primarily HOMs) and/or MM dispersive waves. Given that the spectrum measurements, including the Ge transmission, suggest ≈ 40% of the launched energy is transferred into the Raman soliton, the autocorrelation implies a spatiotemporal distribution similar to those observed in simulations: the Raman soliton occupies lower-order, more-localized modes than the remnant soliton(s) and dispersive radiation. With 21 nJ launched, the 1740-nm Raman soliton has ≈ 80 fs duration and 11 nJ energy. The beam profiles measured after the GRIN fiber show the initial spatial focusing corresponding to the Raman amplification of the lower-order modes (Fig. 6, middle). In a similar experiment with 95 m of fiber, we place a bandpass filter at 1620 nm and adjust the input energy so that the Raman soliton coincides with the filter (spectrum Fig. 6, top right). We thus image the Raman soliton's beam profile (Fig. 6, bottom right) separately from the entire field's beam profile (adjacent). The result confirms that the Raman soliton contains relatively low-order mode content but remains multimode, having a spatial dimension R g = 10 μm = 1.54R o . We note that, as in simulations, the spatial profile of the Raman soliton depends on the initial mode excitation. In some cases, the Raman soliton appears very nearly single mode (R g = R o ). While further work is required to understand the initial conditions necessary for this spatiotemporal fission outcome, the possibility of deterministically controlling the spatial shape of the Raman soliton is promising since it could allow diffraction-limited beams which may be necessary, for example, in microscopy. In agreement with simulations, the Raman soliton remains MM (its spatial width is larger than the fundamental mode by a factor of 1.54) but is more biased to lowerorder modes. The results all display good qualitative agreement with simulations. All scale bars correspond to 20 μm.

Discussion
Several aspects of MM solitons merit further study, such as whether qualitatively different kinds of MM solitons can exist, how they propagate under conditions with weak mode coupling or mode-dependent loss, and how MM solitons collide or interact with one another. Since most other multimode fibers have larger modal dispersion than GRIN fiber, it will be important to establish if multimode solitons can occur within non-degenerate modes in other types of fiber, such as step-index fiber. Although not considered here, promising theoretical results have been achieved for fibers with high disorder [50,51], suggesting MM fibers should provide a platform for the study of incoherent spatially-multimode solitons. In gas-filled hollow-core photonic crystal fibers [25], where MM nonlinear effects have been observed [24, 25], we can furthermore expect to find MM solitons and soliton interactions. Of particular interest will be the fission of MM solitons in the absence of Raman scattering, and the influence of plasma interactions. Finally, we note the significance of this work towards the development of shortpulse fiber sources and other applications. Multimode fibers allow considerably higher pulse energy and average powers than single-mode fibers and, furthermore, MM solitons may have even higher energy than single-mode solitons of comparable spatiotemporal dimensions. As a result, MM solitons may be of key importance in future high-energy ultrashort-pulsed fiber sources. In addition, due to their spatial complexity they may be useful for optical information processing and telecommunications. On the other hand, their spontaneous formation could become a limiting factor in future spatial-division multiplexing schemes.

Conclusion
We have observed multimode optical solitons in GRIN fiber and studied their dynamics. There is good qualitative and often quantitative agreement between simulations and experiments, allowing us to reach the following key conclusions: • MM solitons are observed over a range of spatiotemporal sizes up to ≈ 3 times the fundamental mode size (containing up to an estimated 8-13 modes).
• MM pulses undergo a spatiotemporal soliton fission process, which results in MM solitons and dispersive waves with different distributions of modes.
• During fission, intrapulse Raman scattering typically causes energy from HOMs to transfer into lower-order modes.
• The most intense MM soliton produced from spatiotemporal fission tends to experience a soliton self-frequency shift. The center wavelength of the pulse in each mode redshifts synchronously, so that the entire MM soliton is redshifted.
Ultimately however, much of the behavior and properties of MM solitons remains mysterious. Furthermore, the MM solitons considered in this work represent only a small fraction -coherent multimode solitons in anomalous dispersion GRIN fiber -of the vast parameter space of nonlinear dynamical objects which are expected to be possible in multimode optical fibers.

Appendix
If the field contains a higher-order mode background (HOMB) and a fundamental mode soliton, we can determine how much energy is actually contained in the fundamental mode compare to the total energy, E p . For the energy in the single-mode soliton, Eq. (3) is obeyed: where T (R o , R g ) is the overlap integral between the measured beam and the fundamental mode, E p is the total energy, spread over all modes, and E p T is the energy in the fundamental mode. Meanwhile the relationship between the total energy and the measured pulse duration is: For a Gaussian field E g with waist R g , the overlap integral with the fundamental mode (E o ) can be easily calculated, Dividing Eq. (5) by T , plugging in the above expression and equating with Eq. (6) gives an expression for the relationship between the spatial size of the whole beam, R g and M FHB , For a MM soliton with an actual average size R g , a HOMB can cause the measured size to be inflated to R g . In this case, the HOMB slightly increases the slope measured from the actual value for the MM soliton, M . Repeating the analysis above with R o replaced by R g and M 1 replaced by M yields: