Nonlinear optics in multimode optical fibers (MMFs) is an emerging research field, as it leads to new possibilities for the control of the spatial, temporal, and spectral properties of ultrashort light pulses [1]. Differently from the case of single-mode fibers, which are limited by their small transverse core size, large-area MMF permits for scaling up by orders of magnitude their energy transport capabilities. As a result, research in MMFs has attracted a growing interest for a variety of technologies, e.g., high-power fiber lasers [1], supercontinuum light sources [2], high-resolution biomedical imaging [3], and micromachining [4]. From a fundamental viewpoint, the high beam intensity that can be reached in MMFs has also led to unveiling different complex nonlinear phenomena [5,6].

The optical power transmission of MMFs is limited by different nonlinear effects, whose relative relevance strongly depends on the time scale of the propagating pulses. As well known, in the continuous-wave (cw) (or quasi-cw) regime, the power transmission limitation of optical fibers is set by Brillouin scattering [7]. When subnanosecond sources are employed, different nonlinear loss effects may become relevant. The electrons promoted from the valence to the conduction band by the absorption of pulses longer than tens of picoseconds interact with the atomic lattice, resulting in the heating of the fiber material. Injecting intense light pulses in optical fibers leads to thermally induced irreversible damage, such as fused silica melting and boiling [8,9]. Whenever the temperature reaches extreme values, typically above thousands of Kelvin, the fiber fuse phenomenon takes place [10]. On the other hand, when intense femtosecond pulses are injected into MMFs, the fiber may break down due to ionization phenomena: these typically occur when the laser intensity is above ${10}^{12}\text{W/}{\mathrm{cm}}^2$ [11]. Under proper conditions, photoionization leads to material modifications and ablation. These phenomena have been widely exploited over the past decades, mainly for their application to glass micromachining, e.g., for fabricating fiber Bragg gratings [12].

Owing to the presence of a bandgap, ionization mechanisms in dielectric materials are triggered by multiphoton absorption (MPA) processes [13,14]. MPA plays a significant role even below the breakdown threshold, and it represents a major drawback for optical beam delivery, because it limits the efficiency of optical elements in the high-power regime. In fused silica, which is the main constituent of commercial MMFs, MPA may lead to nonlinear contributions to both the refractive index and the absorption coefficient, when operating at wavelengths between the visible and the mid-infrared range [15,16]. Because of the high intensities associated with the occurrence of MPA, its presence is generally neglected in a telecom context. Indeed, many photons need to be simultaneously absorbed to fill the pure silica bandgap ($>10\mathrm{eV}$) at telecom frequencies: up to six-photon absorption has been observed at $\lambda =790\mathrm{nm}$ [17].

However, the presence of material defects, whose absorption band peaks at a few electronvolts, lowers the number of photons needed for the observation of MPA. For undoped core silica MMFs, the dominant defects, i.e., those with the highest oscillation strength, are the so-called non-bridging oxygen hole centers (NBOHCs) [18]. Recently, it has been shown that pulses at $\lambda =1030\mathrm{nm}$ lead to the simultaneous absorption of 5 photons by NBOHCs. This is followed by the generation of a visible photoluminescence (PL) at the self-focusing point, where the intensity reaches its peak values [19].

In the particular case of graded-index fibers (GIFs), the germanium doping used to shape the core refractive index widens the range of defects, whose oscillation strength is comparable with that of NBOHCs. These defects are the Ge-related oxygen deficiency centers (Ge-ODCs), whose main characteristic is their blue-violet luminescence, and the paramagnetic Ge(1) and Ge(2) centers, which do not show any PL emission [18]. The PL of GIF defects permits directly visualizing the spatial self-imaging (SSI) phenomenon, the latter being the periodic replication of the electromagnetic field upon its propagation [20,21]. SSI is particularly studied in GIFs, since its period is remarkably short (a few hundreds micrometers), owing to the parabolic refractive index shape, which leads to equally spaced propagation constants for nondegenerate modes [22]. As a result, a laser beam propagating inside a GIF continuously widens and tightens its waist: correspondingly, PL generates a periodic array of light spots [23].

In this work, we investigate ultrashort pulse propagation in MMFs, in a power regime that is close but still below their breakdown value. Here MPA mechanisms turn out to be relevant, but no damages of the fiber material are induced yet. This permits us to carry out a detailed experimental characterization of optical nonlinear losses (NLs), which occur both in the normal and in the anomalous dispersion regimes of propagation. We first verify that thermal effects are not involved in determining such losses. Next, by comparing NL measurements of GIFs and step-index fibers (SIFs), we point out that SSI has a key role in enhancing those losses. We investigate the origin of NLs by comparing the wavelength dependence of the nonlinear fiber transmission properties with that of defects’ PL intensity. More specifically, we study the fiber power transmission properties when varying the fiber length, as well as the input beam wavelength and pulse width. In the femtosecond regime, we found that MPA, which is responsible for PL, is one of the main contributors to NLs. To reproduce our observations, we introduce an analytical model, where an effective N-photon absorption term is added to the propagation equation. A value of $N\simeq 3$ is obtained, in good agreement with the experimental transmission at $\lambda =1030\mathrm{nm}$. Our results are important, as they may provide a guideline for modeling the power scaling of spatiotemporal ultrashort pulse propagation in MMFs.

The experimental setup to study NL phenomena in MMFs is shown in Fig. 1. It consists of an ultrashort laser system, involving a hybrid optical parametric amplifier (OPA) of white-light continuum (Lightconversion ORPHEUS-F), pumped by a femtosecond Yb-based laser (Lightconversion PHAROS-SP-HP), generating pulses at 10 to 100 kHz repetition rate and $\lambda =1030\mathrm{nm}$, with Gaussian beam shape ($M^2=1.3$). The pulse shape was measured by using an autocorrelator (APE PulseCheck type 2), resulting in a sech temporal shape with pulse width ranging from 40 to 90 fs for the OPA and between 0.174 and 7.9 ps for the main source. The OPA produces pulses at wavelengths between 650 and 940 nm (signal) and above 1160 nm, up to a few micrometers (idler). As shown in Fig. 1, the laser beam was focused by a 50 mm lens into the MMFs, with an input diameter at $1/e^2$ of peak intensity of approximately 30 μm at 1030 and 1550 nm and of 43.7 μm for the signal beam (between 650 and 940 nm). Both standard 50/125 SIFs and 50/125 GIFs of variable length have been employed. The GIF core radius, core refractive index along the axis, and relative core-cladding index differences at $\lambda =1030\mathrm{nm}$ are $r_c=25\mathrm{\mu m}$, $n_0=1.472$, and $\mathrm{\Delta }=0.0103$, respectively. The input tip of the fiber was imaged by a digital microscope (Dinolite-AM3113T), and light scattered out of the fiber was collected by means of a convex lens into a miniature fiber optics VIS-IR spectrometer (Ocean Optics $\mathrm{USB}2000+$) working between 170 and 1100 nm. The inset of Fig. 1 shows a typical microscope image of a GIF input tip when shining its characteristic PL. At the fiber output, a microlens images the out-coming near field and projects it on a VIS-IR camera (Gentec Beamage-4M-IR). Besides monitoring the output transverse intensity distribution, the CCD also helps to optimize the input coupling. By a cascade of flipping mirror, the beam is then focused into an optical spectrum analyzer (OSA) (Yokogawa AQ6370D) and a real-time multiple octave spectrum analyzer (Fastlite Mozza) with wavelength ranges of 600–1700 nm and 1000–5000 nm, respectively. Finally, the input and output average powers were measured by a thermopile power meter (GENTEC XLP12-3S-VP-INT-D0). To check the presence of fiber damages induced by the propagation of too intense beams, we used an optical microscope (Zeiss Axio Scope.A1).

Optical NLs consist of the transmission decrease when high-power pulses are injected into the fiber. This is shown in Fig. 2(a), where we report the input energy ($E_{\mathrm{in}}$) dependence of the output pulse energy ($E_{\mathrm{out}}$) for a 1 m long GIF. The pulse energy is calculated by dividing the measured average power by the laser repetition rate. We investigated femtosecond pulse propagation both in the normal and in the anomalous dispersion regimes, by using the pump at 1030 nm and the idler at 1550 nm, respectively, whose spectra are reported in Fig. 2(b). As Fig. 2(a) shows, nearly the same behavior was found in both regimes. Whenever $E_{\mathrm{in}}<80\mathrm{nJ}$ (which corresponds to 0.065 and $0.176\text{TW}/{\mathrm{cm}}^2$ of intensity at 1030 and 1550 nm, respectively), these measurements show that $E_{\mathrm{in}}$ and $E_{\mathrm{out}}$ are linearly proportional. This regime is characterized by the presence of weak linear losses, which are of the order of a few decibels per kilometer (dB/km) at 1030 nm [7]. Therefore, for short fibers (such as the case in our experiments) losses due to material linear absorption can be neglected. For simplicity, we do not display injection losses due to coupling misalignment, so that 100% transmission characterizes the linear regime. Figure 2(a) shows that, when $E_{\mathrm{in}}>80\mathrm{nJ}$, the output energy no longer scales linearly with input energy, resulting in a transmission drop, the hallmark of MPA-induced NLs.

In Figs. 2(c) and 2(d), we report the output spectra corresponding to the data of Fig. 2(a). As we can see, the two dispersion regimes lead to strongly different behaviors. At $\lambda =1030\mathrm{nm}$, the output spectra broaden, but they keep peaking at the source wavelength [see Fig. 2(c)]. Conversely, Fig. 2(d) shows that at $\lambda =1550\mathrm{nm}$, a more complex dynamics occurs, leading to a significant enhancement of the nonlinear spectral broadening when compared with the case of Fig. 2(c). As a result of fission of the initial multisoliton pulse, Raman solitons are generated, which are strongly affected by MPA [24]. For wavelengths shorter than 1550 nm, anti-Stokes sidebands arise from spatiotemporal multimode soliton oscillations and almost cover the entire visible range.

By comparing Figs. 2(a), 2(c), and 2(d), we found that, in both dispersion regimes, MPA properties are remarkably similar, in spite of the profound differences in spectral and temporal dynamics. This permits us to justify the derivation, in Section 3, of a simple model for describing the effects of MPA on beam propagation, which only takes into account spatial effects.

When operating at pulse intensities of the same order of magnitude as the laser-induced breakdown threshold, one naturally wonders whether irreversible processes may occur in the fiber. We verified that our experiments can be repeated several times without observing any alteration of the fiber properties. However, we investigate the possible role played by thermal effects. These deal with the vibration of the atomic lattice (phonons), whose motion is much slower if compared to the laser pulsewidth. Therefore, thermal phenomena necessarily must depend on the beam average power, so that their effect of the optical transmission can be experimentally investigated by varying the laser repetition rate. We checked that the fiber output power linearly scales with the laser repetition rate at peak powers in both the linear loss and the NL regimes (Fig. 3). This proves that thermal effects are negligible in our experiments. We also monitored the output transmission during the whole experiment, without noticing appreciable variations. On the other hand, we repeated the experiments with longer pulses that are above 1 ps in temporal duration. Differently from the femtosecond regime, at the same values of peak power, we observed that the fiber transmission slowly but significantly drops in time. In particular, we monitored the transmission for several minutes, observing its progressive reduction (not shown), which is a sign of fiber breakdown. The input tip of the fiber was imaged by the optical microscope working in cross-polarizer configuration in order to maximize the scattering from damages (see the inset of Fig. 3). As can be seen, irreversible modifications of the sample were formed close to the input facet of the fiber, in proximity of the self-focusing point.

The results shown so far exclude the role of dispersion and thermal effects on the physical mechanism responsible for the femtosecond NL. Here we show that the spatial evolution of the beam inside the fiber is the key element for activating the nonlinear optical attenuation. For doing so, we compare the nonlinear transmission properties of MMFs with the same core/cladding size, with (GIF) or without (SIF) SSI. As a result of MPA, one observes the upconversion PL of silica defects, which, thanks to its typical violet color, helps track the spatial dynamics of the beam inside the MMFs. This permits us to visualize the different beam evolutions that take place inside the SIF and the GIF. Optical beams propagating in an SIF experience a single self-focusing point at the very beginning of the fiber [see the digital microscope image of PL in Figs. 4(a) and 4(c)]. The resulting high intensity leads to MPA, so that most of the beam energy is lost over the first few millimeters of propagation. Conversely, in a GIF the beam diameter periodically oscillates in space, owing to SSI [1]. The associated spatial beam breathing produces several minima for the beam diameter, where MPA may lead to NL. As a result, as shown in Figs. 4(b) and 4(d), PL appears as an array of equally spaced emitting points [23]. The different beam dynamics for the two types of MMF have direct consequence on power transmission measurements. In Fig. 4(e), we report the optical transmission as a function of the input pulse energy, for input pulse durations of either 174 fs or 7.9 ps, with $\lambda =1030\mathrm{nm}$. In order to make a meaningful comparison between the two types of MMF, we normalized to 100% the transmission of each fiber at low input pulse energies. In this way, we filter out all the possible linear contributions to losses, which may arise from material absorption or from the scattering at the self-focusing point close to the input fiber tip. As it can be seen, for both MMFs the transmission curves remain flat when picosecond pulses are used, i.e., when the input peak power is not high enough to trigger MPA. On the other hand, using femtosecond input pulses leads to observing strong NL. Moreover, the threshold for a transmission drop is quite different for the two MMFs. As shown in Fig. 2(e), a drop occurs above 80 nJ of input pulse energy in a GIF, whereas a fourfold threshold enhancement is observed in an SIF. This means that the transmission in the presence of NLs is higher in SIFs than in GIFs, whereas the latter are well known to better confine and propagate optical beams in the absence of MPA. For this reason, in order to highlight the role of NLs, only experiments carried out with GIFs will be discussed in the following subsection.

With the aim of investigating the physical origin of MPA, we studied how both PL and fiber transmission vary with the source wavelength. As depicted in the setup of Fig. 1, we collected the PL signal into a VIS-IR spectrometer. In Figs. 5(a) and 5(b) we show measured scattered spectra and input spectra, respectively, for a laser peak power of 2.5 MW, and wavelengths varying between 680 and 900 nm. Besides the spectral broadening of the pump, observed spectra display the characteristic NBOHC and the Ge-ODC PL peaks, occurring at 650 and 400 nm, respectively. Following the same method as in Ref. [19], we suppose a power-law linking the Ge-ODC PL intensity ($I_{\mathrm{PL}}$) and the input peak power ($P_p$): $I_{\mathrm{PL}}\propto P_p^{N_{\mathrm{PL}}}$, where $N_{\mathrm{PL}}$ is the average number of photons involved in the MPA process exciting the defects PL. Thus $N_{\mathrm{PL}}$ can be obtained by evaluating the slope of the linear dependence from $P_p$ of $I_{\mathrm{PL}}$ (calculated as the integral of the corresponding peak in the spectrum) when displayed in a log–log plot. As reported in Fig. 5(c), the value of $N_{\mathrm{PL}}$ varies with the source wavelength. Specifically, about 3 photons are simultaneously absorbed at $\lambda =750\mathrm{nm}$, since the Ge-ODC absorption band has a maximum at about 250 nm [18]. We do not report here the corresponding analysis for NBOHC, since its contribution to PL is so small that the resulting linear fit is too noisy to provide a reliable estimation of $N_{\mathrm{PL}}$. Furthermore, the NBHOC PL signal partially overlaps with the source at small wavelengths, as can be seen in Fig. 5(a).

Next, we compared the analysis of the wavelength dependence of side-scattered spectra with that of fiber transmission measurements. As shown by Fig. 5(d), NL was observed at all source wavelengths. Nonetheless, differently from the PL related curves of Fig. 5(c), the transmission curves are not sorted with respect to the input wavelength. This indicates that Ge-ODC and NBOHC absorption is not the only mechanism that is responsible for NLs. For a clearer comparison, in Fig. 5(e) we simultaneously plot the wavelength dependence of both $N_{\mathrm{PL}}$ and the transmission at $P_p=4\mathrm{MW}$. As can be seen, while $N_{\mathrm{PL}}$ monotonically increases with the source wavelength, the transmission curve exhibits a complex behavior. Specifically, two transmission maxima were obtained around 750 and 840 nm, respectively. The first maximum corresponds to 3 times the Ge-OCD and NBOHC absorption band peaks, whereas the second maximum can be ascribed to other types of defects that do not contribute to PL, such as the Ge(1) and Ge(2) centers, whose associated oscillation strengths are comparable to those of luminescent defects [18].

In order to track the spatial evolution of PL along the beam propagation, we placed the spectrometer and the lens on a translation stage, which allows for shifting them solidly. In Fig. 6(a) we show the longitudinal evolution of the side-scattered light spectrum, over the first centimeter of a GIF, for 2.5 MW of input peak power. Besides the pump at 1030 nm, three additional main peaks appear: the third harmonic (TH) at 343 nm and two PL peaks at 400 and 650 nm, respectively. The resulting collected spectra were postprocessed in order to determine the evolution of the spectral energy associated with each peak. Figure 6(b) shows that the energy of PL peaks significantly drops over less than 1 cm of fiber. On the other hand, our measurements (not shown here) reveal that the TH (and the pump) energy is only damped over several centimeters (or several tens of centimeters). The observed decrease of the PL is associated with a drop of the peak intensity of the beam. However, the drop of PL intensity does not necessarily prove a loss of the beam energy. For example, an intensity decrease may be due to a temporal broadening of the pulses, to modal walk-off, or to an increase of the effective area of the beam. However, the observed rapid damping of light side-scattered by the pump indicates that, even if present, these contributions are negligible in our case [i.e., for a propagation distance of the order of 1 cm, as in Fig. 2(a)].

In order to verify that indeed MPA processes are responsible for the observed PL signal drop of Fig. 6(b), we performed a cut-back experiment by progressively reducing the GIF length from 10 cm down to 1.5 cm, while keeping the input coupling conditions unchanged. The measured transmission as a function of $P_p$ is reported in Fig. 6(c). $P_p$ was varied by changing the laser compressor ratio. Namely, thanks to a feedback control of the laser system, the energy of each pulse was kept a constant, while the pulse duration was varied between 174 fs and 7.9 ps. The advantage of this method is that optical elements that may introduce small misalignments, such as variable attenuators, are not needed. As previously observed with reference to Fig. 2(a), two distinct power transmission regimes are revealed. For powers below the threshold value of $P_p=0.5\mathrm{MW}$, the fiber transmission remains nearly constant for all fiber lengths; whereas at higher powers, the transmission starts decreasing. It is interesting to point out that such a power threshold coincides with the value that leads to the appearance of PL [see Fig. 6(d)]. This observation suggests that MPA leading to PL is one of the main mechanisms that are responsible for NLs. However, other nonlinear effects, such as the aforementioned nonluminescent defect excitation, may give a contribution to NLs. In the next section, we propose a model based on an effective N-photon absorption, which is able to quantitatively reproduce the experimental transmission drop.

Let us consider the beam spatial dynamics under the slowly varying envelope approximation, treating time as a parameter, analogously to the variational approach proposed by Karlsson et. al. [25]. In order to take into account the experimental NL, we added to the index grading and the Kerr nonlinearity, an effective $N$-photon absorption term into the GIF permittivity: $$\epsilon \simeq n_0^2+2n_0{n}_2{|E|}^2-g{n}_0^2{r}^2+i{n}_0\frac{{\alpha }_N}{k_0}{|E|}^{2N-2}.$$(1)Here, $n_2$ is the nonlinear refractive index, $E$ is the envelope of the propagating electric field in the fiber (measured in $\sqrt{\mathrm{W}/\mathrm{m}}$), $g$ is the index grading parameter defined as $g=2n_0\mathrm{\Delta }/r_c^2$, $r$ is the independent radial coordinate, $k_0=2\pi n_0/\lambda$ is the propagation constant in the core, and ${\alpha }_N$ is the N-photon absorption coefficient. We highlight that $N$ represents the average number of photons involved in the nonlinear absorption process, which sums all possible contributions to energy losses, including side-scattering. As in Ref. [25], we impose a Gaussian shape of the beam, so that the electric field can be written as $E(z,r)=A(z)e^{-\sigma (z)r^2}$. However, in a nonconservative system, the Lagrangian equations of the variation approach do not hold, and the beam power defined as $$P=2\pi \int {|E|}^{2}r\mathrm{d}r$$(2)is no longer an integral of motion. Here, we consider an alternative method to the Lagrangian equations: we impose that only the $z$-dependent amplitude $A$ is affected by absorption, so that $\sigma$ is independent of ${\alpha }_N$. In this way, we can recover the evolution of $\sigma$ from the lossless problem. With this consideration, the beam power evolution equation reads as (see Appendix A for the step-by-step derivation) $$\frac{\mathrm{d}P}{\mathrm{d}z}=-\frac{{\alpha }_N}{2N-1}\frac{P^N}{{(\pi a^2)}^{N-1}},$$(3)where $a(z)$ is the beam $1/e^2$ radius that oscillates along $z$ due to the SSI, according to the formula $$a(z)=a_0\sqrt{{\cos }^2(\sqrt{g}z)+C{\sin }^2(\sqrt{g}z)}.$$(4)The parameter $C$ contains the dependence of the beam size upon input peak power, and it is defined as $$C=\frac{1-p}{{k_0}^2{a}_0^{4}g},$$(5)with $p=n_2{k_0}^2{P}_p/2\pi n_0$ and $P_p=\pi A_0^2{a}_0^2$. The fiber breakdown condition is given by imposing $C=0$, which corresponds to the peak power threshold $P_p^{\mathrm{t}\mathrm{h}\mathrm{r}.}=4.26\mathrm{MW}$ (using the experimental parameters at $\lambda =1030\mathrm{nm}$ and $n_2=2.7\times {10}^{-20}{\mathrm{m}}^2/\mathrm{W}$). While the SSI period is power independent, the beam size oscillation amplitude depends on $P_p$ [as depicted in Fig. 7(a)] [23]. In Fig. 7(b), we report the oscillating behavior of the beam size along the propagation distance, for the power $P_p=2\mathrm{MW}$, right below the breakdown threshold.

In an MPA process, the energy loss strongly depends on the beam waist through its intensity. At points of maximum intensity, the MPA contribution becomes most relevant. Even though the absorption coefficient is small, the intensity locally reaches such high values that the NL term in Eq. (3) is no longer negligible. This results in a series of stepwise drops of the transmission curve at points of SSI-induced beam focusing as shown in Fig. 7(c).

One may object that $a$ should depend on $P$, instead of its initial value $P_p$, since the presence of absorption modifies the total energy of the beam along its propagation. Nevertheless, we numerically verified that considering $a$ dependent on $P_p$ instead of $P$ does not introduce a significant overestimation of NL. As a matter of fact, ${\alpha }_N$ keeps the same order of magnitude when using the experimental parameters, regardless of the fact that $a$ is assumed to depend on $P$ or $P_p$. This observation provides a useful way for simplifying the fitting of the experimental data, which is the topic of the next section.

As a last result, we fitted the experimental data of the cutback measurements in Fig. 6(c) with our analytical model (see Appendix B for details). As shown in Fig. 8, we obtained a good quantitative agreement: the fitting function interpolates well experimental data for all values of $P_p$ and $z$. In Fig. 8, we juxtapose 2D and 3D plots of the transmission in order to appreciate the goodness of the fit. The extrapolated parameters are $N=3.008$ and ${\alpha }_N=2.415\times {10}^{-33}$ (in SI units). As already foreseen in the comparison between the transmission curves and $N_{\mathrm{PL}}$, a value of $N\simeq 3$ indicates that MPA that excites PL is not solely responsible for the observed NL. Many other factors contribute to the value of $N$, such as Ge(1) and Ge(2) defect absorption and generation of the pump TH, which is not guided [and therefore it is laterally scattered from the fiber as shown in Figs. 5(a) and 5(c)]. As a general remark, our model is rather simple, since it considers a monochromatic wave, in spite of the broad spectra of Fig. 2(c). Furthermore, we are not considering effects such as multiphoton ionization [9], color-center generation [26], tunneling photoionization [14], and plasma formation [17], which have been shown to play a role when MW peak power femtosecond lasers are employed. Models of femtosecond laser absoption in silica, which take into account filamentary propagation, ablation, and refractive index micro-modification, can be found in the literature [8,27,28]. Nevertheless, we believe that our model achieves its purpose, as it finds a remarkably good agreement with the experimental data. Specifically, we could reveal that an effective three-photon absorption term is useful for describing the beam dynamics in a variational-like approach when working close to the fiber breakdown threshold.

We experimentally characterized nonlinear optical losses in MMFs when operating slightly below the breakdown power threshold. In this regime, MPA mechanisms play a significant role, and may strongly affect the dynamics of the propagating beam. For example, in the anomalous dispersion regime, it has been shown that MPA clamps the beam output energy, thus suppressing the Raman soliton self-frequency shift [24]. In this work, based on the observed close similarity between nonlinear power transmission properties of MMFs in both the normal and the anomalous dispersion regimes, we carried out a detailed investigation of the spatial beam dynamics. Specifically, we revealed that SSI has a key role in enhancing NL. Therefore, MPA effects are remarkably higher in GIFs than in SIFs. By comparing side-scattering and transmission measurements, we could infer that the observed NL cannot be entirely ascribed to the MPA, which is associated with defects’ PL. Although all of the different sources of NL could not be directly identified, we demonstrated that the observed nonlinear transmission drop could be quantitatively well reproduced by introducing an effective N-photon absorption term in the propagation equation. Namely, we found that $\mathit{N}=3$ at $\lambda =1030\mathrm{nm}$ provides a good agreement with experiments carried out at several input peak powers and fiber lengths. Our results highlight the intrinsic limitation of VIS-IR light propagation in MMFs due to MPA, which becomes highly relevant when propagating femtosecond pulses of MW peak powers, i.e., close to, but still below, the fiber breakdown threshold. In this sense, our results will be of significant interest for different applications based on MMFs, such as micro-machining, medical imaging, beam delivery, and fiber lasers for spatiotemporal mode locking. In all of these emerging technologies, MPA represents a major obstacle for their power upscaling.

Acknowledgment. We acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program, and the Italian Ministry of University and Research.