Nonlinear beam self-imaging and self-focusing dynamics in a GRIN multimode optical fiber: theory and experiments

Beam self-imaging in nonlinear graded-index multimode optical fibers is of interest for many applications, such as implementing a fast saturable absorber mechanism in fiber lasers via multimode interference. We obtain an exact solution for the nonlinear evolution of first and second order moments of a laser beam carried by a graded-index multimode fiber, predicting that the spatial self-imaging period does not vary with power. Whereas the amplitude of the oscillation of the beam width is power-dependent. We have experimentally studied the longitudinal evolution of beam self-imaging by means of femtosecond laser pulse propagation in both the anomalous and the normal dispersion regime of a standard telecom graded-index multimode optical fiber. Light scattering out of the fiber core via visible fluorescence emission and harmonic wave generation permits us to directly confirm that the self-imaging period is invariant with power. Spatial shift and splitting of the self-imaging process under the action of self-focusing are also emphasized.


INTRODUCTION
Nonlinear multimode optical fibers (MMFs) are an emerging research field, as they permit new ways for the control of spatial, temporal and spectral properties of ultrashort pulses of light [1]. As a result, nonlinear MMFs are of interest for a diversity of optical technologies, e.g., for scaling-up the power of fiber lasers and supercontinuum light sources, for high-resolution biomedical imaging, micromachining, and high-power beam delivery, to name a few. MMFs also provide a simply accessible testbed for the study of complex physical phenomena [2]. Among MMFs, graded-index fibers are of particular interest for nonlinear optics studies, since the reduced modal dispersion leads to relatively long interaction lengths among the fiber modes, even for pulses in the femtosecond regime [3].
Spatial self-imaging (SSI), first observed by Talbot in 1836 [4], is a peculiar property of beam propagation in GRIN fibers. SSI leads to periodic oscillations along the fiber of the beam width and intensity. In combination with the Kerr effect, SSI leads to a longitudinal periodic modulation of the refractive index of the fiber core, akin to a dynamic long-period grating. As a result, as recently discussed in a review paper by Agrawal [5], SSI leads to many recently discovered nonlinear effects in GRIN MMFs, such as dispersive wave series emission from multimode femtosecond solitons [6,7], geometric parametric instability of Continuous Wave (CW) beams in the normal dispersion regime [8], and spatial beam self-cleaning phenomena [9,10]. In addition, SSI has been recently widely exploited for the mode-locking of fiber lasers, by exploiting multimode interference (MMI) resulting from SSI in a short piece of GRIN MMF with precisely controlled length, sandwiched between two singlemode fibers.
Because of SSI, any input field profile is periodically reproduced at equally spaced points z s along the GRIN fiber, such that where β n and β 1 are the propagation constants of higher-order modes with index n and of the fundamental mode, respectively, and m n is an integer. As a result, any input beam shape is reproduced at the fiber output, whenever its length is exactly equal to an integer multiple of z s . However, unless the shape of the input beam exactly matches the fundamental mode of the fiber, the input beam is spread in several propagating modes. As a result, the periodic reconstitution of the initial beam shape is obtained via the superposition of modes with different propagation constants, which affects the output beam divergence. The presence of self-imaging in the MMF was indirectly put into evidence by Zhu et al. [11] by tuning the wavelength of the input signal, and measuring the transmission spectrum of the MMI structure. Nazemosadat and Mafi proposed to use intensity-dependent differential phase shift among transverse modes, to obtain nonlinear MMI in a short length of GRIN MMF [12]. According to their description, in a nonlinear MMF the self-imaging period should change to z l where (β n (I) − β 1 (I))z I = πm n . ( In a realistic situation, an input laser beam coupled to a GRIN MMF excites several modes, each of them carrying a different amount of intensity along the fiber. Therefore the simple description of Eq. (2) should be substituted by the solution of nonlinear coupled mode equations including self-and cross-phase modulation, as well as four-wave mixing terms [13]. Equivalently, one may use a two-dimensional nonlinear Schrödinger equation (2D-NLSE) [9].
By sandwiching the MMF between two singlemode fibers (SMS structure), it has been experimentally observed that nonlinear MMI acts as a fast saturable absorber (SA) mechanism which permits to obtain mode-locking in high pulse-energy fiber lasers. Specifically, low power signals are strongly attenuated when propagating through the SMS structure, while high intensity pulses experience a relatively high transmission. Nonlinear MMI-based SAs permit to operate at much higher pulse energies and peak powers than other SA mechanisms, thanks to their high damage threshold, low cost, simple structure, and mechanical robustness. Nonlinear MMI in a GRIN MMF also permits spectral filtering, a necessary property for fiber laser mode-locking in the normal dispersion regime.
Several different variants of the technique of fiber laser mode-locking based on nonlinear MMI have been demonstrated in recent years [14][15][16][17][18][19][20][21][22]. In order to test the nonlinear transmission of the SMS structure, the MMF was stretched, and the SMS transmission was measured as a function of the stretching length, for different input intensities [18]. It was found that, at relatively high intensities, the self-imaging induced beam oscillations occur with much smaller amplitude, and with average (i.e., over one self-imaging period) transmission that is significantly increased ( 75%) with respect to the low intensity case ( 68%). This suggests that nonlinear mode coupling and the resulting beam reshaping could be a mode-locking mechanism, rather than the nonlinear variation of the self-imaging period. Indeed, experiments show that the excitation of high-order modes (HOMs) in the GRIN fiber is a necessary condition to obtain saturable absorber action [15,16]. A study of the impact of the GRIN core diameter (ranging from 20 to 62.5 µm) revealed that larger diameters (hence larger proportions of excited HOMs) increase the output laser power. Another mode-locking mechanism could result from the longitudinal translation of the self-imaging process, because of self-focusing effects in the GRIN MMF.
Already back in 1992, Karlsson et al. have theoretically studied, by means of the variational method, the dynamics of self-imaging in a GRIN MMF under the combined action of diffraction, nonlinearity, and parabolic index profile [23]. They derived an approximate analytical solution for a multimode beam, which permit to evaluate the evolution of the beam width and its longitudinal phase delay. In this work, we derive an exact solution for the nonlinear evolution of first and second order moments of a laser beam carried by a GRIN multimode fiber. Our theoretical analysis does not require, unlike methods based on the variational approach, that the beam maintains a specific shape (e.g., gaussian) throughout its propagation in the MMF. This permits us to show, in full generality, and in contrast with current understanding of nonlinear MMI based SAs, that the SSI period does not vary with power. Whereas the amplitude of the beam width oscillations is power-dependent, which is still sufficient to enable a power-dependent transmission for the SA. We also derive the power threshold for catastrophic self-focusing in the GRIN MMF.
Although SSI is a well-known process, it remains difficult to directly prove its existence in MMFs. Here we show that it is possible to work around this problem, by experimentally recording the local nonlinear parametric conversion obtained in a non-collinear geometry, i.e., by using Cherenkov phase matching [24]. As we shall see, the periodic emission of a second harmonic wave, obtained at the core-cladding interfaces, accompanied by multi-photon absorption leading to wideband fluorescence in the blue spectral domain, provides a clear evidence of the periodic evolution of the power density in the MMF. We could highlight these processes, we believe for the first time, by using femtosecond laser pulse excitation both in the anomalous and in the normal dispersion regime. This permits us to obtain sufficiently high peak powers, in the multi-MW range, at the points of minimum beam waist of the SSI process.

Model
We consider beam propagation in a multimode optical fiber with a parabolic index profile. The beam dynamics is assumed to be described by a 2D-NLSE with an instantaneous Kerr nonlinearity of the form where k 0 = ω 0 n co /c is the wavenumber, ρ is the core radius, ∆ = (n 2 co − n 2 cl )/2n 2 co is the relative index difference, and n co (n cl ) is the refractive index of the fiber core (cladding), respectively. Here we neglect dispersion and Raman scattering and assume that the index profile has an infinite extent (i.e. it does not truncate when the cladding index is reached).

Moments of the beam
Now, we show that it is possible to obtain an exact, analytical solution for the moments of the field, whose evolution is described by the 2D-NLSE (3). In particular, we derive a closed system of equations for the various moments. The moments may conveniently be obtained by considering them as expectation values of observables and using operator methods from quantum mechanics. With this machinery the evolution equations for the various moments can be found in a few lines of calculation using commutator algebra. We introduce a linear Hamiltonian operator where β = k 0 ∆/ρ 2 , γ = k 0 n 2 /n co ,x andŷ are transverse position operators, and the nonlinear termÎ = |A(x, y)| 2 is treated as a potential term that should be determined self-consistently. Here we have also made use of the momentum operatorsp x = −i∂/∂x andp y = −i∂/∂y. With this pseudo-Hamiltonian we can write Eq. (3) as Following the procedure described in the Appendix A, we find the sinusoidal solutions for the first-order moments where the oscillation period depends on the index parameters, and a x,y , b x,y are integration constants that depend on the initial conditions. Moreover, the radial moment x 2 + y 2 = x 2 + y 2 , associated with the rms-width for a centered beam, satisfies the closed equation where the Hamiltonian invariant H 0 is defined in Eq. (27). This equation has a solution when H 0 > 0 that describes sinusoidal oscillations around a constant average value, viz.
where a r and b r are integration constants that determine the amplitude and phase of the oscillations. These can be found from the initial beam profile entering the fiber. Whenever the input beam is not in the process of changing, i.e.
[d x 2 + y 2 /dz] 0 = 0, one obtains a r = x 2 + y 2 0 − H 0 /4β and b r = 0 . The oscillation amplitude is consequently a function of the input beam width, while the SSI period is determined by the fiber parameters only, and it remains equal to its linear value z s = πρ/ √ 2∆. It is also seen that the oscillation period of the first-order moments is twice that of the rms-width.
The first integral of Eq. (7) represents a conservation law, M 2 − γk 0 x 2 + y 2 I = const., for the beam quality factor [25]. Furthermore, assuming b r = 0, we may rewrite Eq. (8) by using the beam compression parameter (C-parameter) that was introduced by Karlsson et al. [23] x 2 + y 2 (z) = This parameter measures the relative importance of diffractive broadening and nonlinear self-focusing. The radial moment is seen to be stationary for C = 1, while it becomes negative for C < 0 which corresponds to the condition for collapse. The critical collapse power is reduced in a GRIN fiber with respect to a homogeneous bulk material, because of the guiding refractive index profile. In the latter case the Hamiltonian must be negative for collapse to occur for a beam without an initial phase front curvature [26].

Vortex beam
A numerical example of a propagating beam with a non-trivial beam profile is shown in Fig. 1, together with the corresponding evolution of the associated moment. The normalized simulation parameters are k 0 = 1.3, ∆ = 1.7, ρ = 0.8, n co = 2.0 and n 2 = 1.5, and the initial field is given by the function A = A 0 sech(x/a x )sech(y/a y ) sin(x)e i2tan −1 (y/x) with coefficients a x = 0.81, a y = 1.01 and A 0 = 0.56 that has been shifted to the position (x 0 , y 0 ) = (1.29, 2.47). The simulated moment (red dots) is found to be in excellent agreement with the analytical expressions that is shown as a blue curve in Fig. 1 (any discrepancy is only due to numerical inaccuracy).

Gaussian beam
It is also evident from Eq. (8) that the constant average value for the rms width depends on the Hamiltonian invariant that is determined by the initial conditions. Assuming a Gaussian initial beam profile of the form The (1/e) beam width is found by calculating the radial moment x 2 + y 2 / P 0 = (x 2 0 + y 2 0 )/2, which reduces to x 0 = y 0 for a symmetric beam. It is seen that the magnitude of the Hamiltonian is reduced by the nonlinear term for a focusing nonlinearity. This implies that the average value of the radial beam width will shrink as the power is increased. In Fig. 2 we illustrate the power-induced variation of the average beam width, i.e. w = H 0 /4β P 0 , for different values of the input width x 0 . As can be seen, the average beam width decreases as the power P 0 grows larger: the narrower the input beam size x 0 , the faster the nonlinear decrease of the beam width with power.
When using fiber parameters from the numerical simulations reported in Ref. [9], and assuming a fixed input beam diameter of 40 µm (FWHM) so that x 0 = y 0 = 24.0 µm, one obtains that associated power-induced variation of the average width for the radial moment w occurs on a GW scale (and w shrinks down to zero for P 0 = 2.5 GW). This means that the nonlinear reduction of the Hamiltonian is not the root cause of spatial beam cleaning. Note that, in this case, the average radial beam width approaches the value w = 17.0 µm (c.f. core radius ρ = 26 µm) in the low power limit. This can be compared with the radial width of the fundamental linear eigenmode solution which has w f = x 2 + y 2 / P 0 = 1/ √ α ≈ 5.76 µm or 9.59 µm FWHM (see corresponding orange curve on Fig. 2). Note that, for an input beam perfectly matching the fundamental mode shape, the average beam width shrinks until it reaches zero when P 0 ≈ 16.55 MW. However, the beam will not be stationary in the GRIN profile, and oscillations of the beam width will cause it to collapse at a lower critical power. The moment theory is however exact under the assumption of a non-truncated index profile, which suggests that a cleaned beam should still have an overall rms-width, although perhaps with a skewness, that is in agreement with the theory. To better highlight the nonlinear dependence of the solution for the second-order radial moment Eq. (8), we plot in Fig. 3 the input beam power P 0 dependence of the beam width oscillation amplitude A r = 4βa r / H 0 , for different values of the the input radial width w. As can be seen, if the initial width w = w f , the oscillation amplitude is equal to zero in the linear limit, and it increases with power until it reaches unity for P 0 8.3M W . This means that the beam width shrinks to zero, which corresponds to a beam collapse condition. In fact, the C-parameter shows that, for a symmetric Gaussian with flat initial phase, the critical collapse power P c = 2π/(γk 0 ) = λ 2 /(2πn co n 2 ) is independent of the beam width. For a beam with a nonvanishing initial phase-front curvature, our approach also permit to obtain the condition for critical collapse to occur. It is interesting to point out that the collapse power depends on the incidence angle of the beam into the GRIN MMF. Details about the analytical description of SSI beam evolution in the most general case is provided in the Appendix. Fig. 3 also shows that, in the low-power limit, the oscillation amplitude is positive (negative) when the input beam width is larger (smaller) than the width of the fundamental mode w f . A positive (negative) value of A r means that the beam width is decreasing (increasing) along the fiber with respect to the input value, until it experiences a maximum beam compression (widening), before returning back to the input width after one period. At sufficiently high powers, self-focusing leads to a narrowing beam width in all cases.
The generality of the moments method also permits to analytically study the SSI dynamics for beams with non-Gaussian transverse profile. We have thus generalized the analytcal description to the case of a super-Gaussian initial beam profile: corresponding results are reported in the the Appendix.

EXPERIMENTS
In order to confirm the theoretically predicted invariance of the SSI spatial period z s with respect to the input power, we carried out an experimental study of the dynamics of SSI in the nonlinear regime of pulse propagation in a GRIN MMF. By using femtosecond input pulses with peak powers up to the MW power range, we could directly visualize the high intensity points inside the MMF, thanks to the associated side scattering of visible higher-harmonics and fluorescence light.

Anomalous dispersion
We performed two sets of measurements, with two different laser sources. First, we used an ultra-short femtosecond laser system, involving a hybrid optical parametric amplifier (OPA) of white-light continuum, pumped by a femtosecond Yb-based laser, generating 120 fs pulses at 1550 nm, with 25 kHz repetition rate. The input laser beam was focused by a 30 mm focal lens, corresponding to an input beam diameter (1/e 2 ) of 18 µm, into a 5 cm long multimode standard 50/125 GRIN fiber, with relative index difference ∆ = 0.0102.   Figure 5, the dependence of SSI as a function of the input peak power of the injected pulses. Here the input peak power that is coupled into the MMF is increased from 2.3 MW up to 7 MW. Figure 4 shows, in the top panel, the picture taken by a digital microscope. The second panel from the top shows that blue fuorescence is periodically scattered from the side of the MMF. The bottom two panels in Figure 4 show that green light is also scattered by the defects of the fiber cladding.
As can be seen in Figure 5, the SSI period remains unchanged, and remarkably close to the theoretically predicted value (i.e., z s = πρ/ √ 2∆ ≈ 550µm). The input pulses generate high-order multimode solitons in the fiber. These solitons undergo, after a propagation distance of about 10 cm, fission into several fundamental solitons under the action of Raman soliton self-frequency shift and higher-order dispersion. Here we limit ourselves to consider the regime of multisoliton propagation over the first few cm of GRIN MMF, that is before that soliton fission takes place. The dynamics of multimode Raman solitons generated by the fission process will be discussed in a separate publication. However, it is important to point out that the observed processes of multiphoton absorption-induced side-scattering of visible fluorescence lead to significant nonlinear losses, which induce intensity clamping at the output of the GRIN MMF.
The invariance of the SSI period with power is well illustrated by Figure 6(a), that shows the input power dependence of the spatial frequency intensity spectrum of the scattered fluorescence. These results confirm that the SSI period remains a constant (within experimental measurement errors) even at the highest input peak powers.
We also measured the power-dependence of the beam size at the points of maximum beam compression, corresponding to the bright spots in Figure 5. Figure 6(b) reports the dimension of the beam in the transverse section of the fiber: as can be seen the beam size does not show a significant dependence on input power. This confirms the hypothesis of Section that the ansatz for the transverse mode profile is maintained along the nonlinear beam propagation. On the other hand, Figure 6(c) shows that the beam size in the axial dimension exhibits a power-dependence.
Specifically, Figure 6(c) shows that, for input peak powers below (above) 4 MW, the beam dimension in the axial direction increases (decreases) along the fiber. The longitudinal dependence of the beam size in the axial direction could be related with the presence of the significant nonlinear losses, which occur over the first few centimeters of the fiber. To put into evidence the nonlinear transmission properties of the GRIN MMF subject to input peak powers just below the critical value for collapse, we reported in Figure 7 the values of the fraction of input coupled energy that emerges from different lengths of the fiber, as a function of the input peak power. As can be seen, from 1 cm of fiber the transmission drops below 50% at the highest peak powers close to 5 MW. Whereas for lengths above 7 cm, the transmission drops below 20% for powers approaching 3 MW. Note that the different curves for fiber lengths above 7 cm tend to overlap, indicating that most of the nonlinear loss occurs over the first few centimeters of the fiber. We also carried out a series of measurements in the normal dispersion regime of the fiber, by using a fiber laser source at 1030 nm, generating 250 fs pulses with a 30 kHz repetition rate. Here we injected the pulses in a relatively short, 1 cm section of a 50/125 µm GRIN MMF. We obtained at the nodes (that is, at the points of minimum beam waist and maximum peak intensity) of SSI a sufficient intensity to trigger, in addition to multi-photon fluorescence, also non-collinear (i.e., using Cherenkov phase matching [24]) second-harmonic generation (SHG), both of which are scattered outside the fiber cladding (see Figure 8). The SHG is due to the presence of a weak quadratic nonlinearity because of the Ge doping ions. Beyond the observation of the self-imaging periodicity by means of non-collinear frequency conversions, we investigated possible longitudinal distortions introduced by a self-focusing regime (see Figure 9). By increasing the input beam power (up to 5.2 MW), and keeping the on-axis excitation, we initiated a self-focusing propagation regime, starting from the first node of the self-imaging process.
This spatial trapping, which resembles a Townes soliton [27], propagates over hundreds of micrometers before recovering its diffracting nature under the effect of nonlinear losses introduced by frequency conversions. Because of its intensity-dependent nature, this extreme event can drastically modulate the initial mode beating, by introducing both a pulse breaking process and a shift of the self-imaging periodicity. As a result, the intense beam of Figure 9(b) is transformed in a double peak of intensity, both oscillating with the same initial period of the self-imaging process.
As shown by Figure 10, light scattering outside the fiber, at 1030 nm, gives again a clear visualization of the presence of SSI, which leads to sharp intensity peaks along the periodic evolution of the field in the GRIN MMF. The measured periodicity of the light intensity in the fiber matches well the value obtained with experiments at 1550 nm, as well as the theoretical value.
We also investigated the case when the input beam is injected at a small angle with respect to the axis of the fiber. In this case, as shown by the resulting series of bright spots of Figure 10, the beam traces out a zig-zag trajectory, as it appears to be reflected by the core-cladding index boundary.
It is also clearly visible that the transverse dimension of the spot does not extend over all the transverse fiber core section, but is either localized on the central part or on a side. Thus, the longitudinal modulation imprinted on the propagating light inside the fiber is directly dependent of the initial coupling conditions, i.e. of the combination of initially excited modes which can self-modulate because of the mode beating and the Kerr effect, and initiate a spatial exchange of energy between them. This mechanism is at the origin of the spatial self-cleaning effect, which has been largely reported in the literature [9,10]. Additionally, this experiment demonstrates that the transverse position of the nodes of the periodic beating can be transversely tuned by means of the initial coupling conditions. This brings an additional degree of freedom to realize a nonlinear saturable absorption mechanism, eventually favouring the emergence of a high-order mode, in order to benefit of the higher dispersion regime when light is coupled back in the multimode fiber. In this sense, temporal mode-locking on a high-order transverse mode should be possible. CONCLUSION We studied the dynamics of beam-self imaging in nonlinear GRIN multimode optical fibers. We obtained an exact solution for the first and second order moments of a laser beam, describing both the period and the amplitude of the beam width oscillations along the fiber. The theory also permits to analytically predict the critical power for beam critical self-focusing, or collapse. We experimentally studied the longitudinal evolution of beam self-imaging by means of femtosecond laser pulse propagation in both the anomalous and the normal dispersion regime of a standard GRIN fiber. By observing light scattering out of the fiber core via visible fluorescence emission and harmonic wave generation, we could directly confirm that the invariance of the self-imaging period up to values close to beam collapse. These findings are of interest for applications involving fiber lasers mode-locked via multimode interference, and to all-optical beam processing with multimode fibers.
where i, j denotes combinations of the x, y operators, and the commutator is defined by [f ,ĝ] =fĝ −ĝf . We assume that = 1 and that the commutators obey the standard Lie algebra To derive the moment equations we make use of the fact that the integrals can be expressed as expectation values of operators q = AqA * ds, where ds = dxdy. As a result, the evolution equation for each moment are obtained from In particular we have that and analogous for y.
where the C-parameter satisfies and is defined by the condition that the minimum rms-width is x 2 + y 2 = C x 2 + y 2 0 . Collapse will therefore occur at some point along the beam trajectory when C ≤ 0.

Supergaussian beam profile
To investigate the influence of the beam shape on the dynamics we also consider a supergaussian beam profile A(x, y, 0) = A 0 exp − 1 2 The power and Hamiltonian invariants for this case are obtain as while the beam width and C-parameter becomes Also here it is seen that the critical collapse power is independent of the beam width for a symmetric beam.