Nonlinear beam self-imaging and self-focusing dynamics in a GRIN multimode optical fiber: theory and experiments
Tobias Hansson, Alessandro Tonello, Tigran Mansuryan, Fabio Mangini, Mario Zitelli, Mario Ferraro, Alioune Niang, Rocco Crescenzi, Stefan Wabnitz, Vincent Couderc
NNonlinear beam self-imaging and self-focusing dynamics in a GRIN multimode opticalfiber: theory and experiments
Tobias Hansson , Alessandro Tonello , Tigran Mansuryan , Fabio Mangini , Mario Zitelli ,Mario Ferraro , Alioune Niang , Rocco Crescenzi , Stefan Wabnitz , , and Vincent Couderc Department of Physics, Chemistry and Biology,Linköping University, SE-581 83 Linköping, Sweden Université de Limoges, XLIM, UMR CNRS 7252,123 Avenue A. Thomas, 87060 Limoges, France Dipartimento di Ingegneria dell’Informazione, University of Brescia, Via Branze 38, 25123 Brescia, Italy Dipartimento di Ingegneria dell’Informazione, Elettronica e Telecomunicazioni,Sapienza University of Rome, Via Eudossiana 18, 00184 Rome, Italy and Institute of Automation and Electrometry SB RAS,1 ac. Koptyug ave., Novosibirsk 630090, Russia
Beam self-imaging in nonlinear graded-index multimode optical fibers is of interest for manyapplications, such as implementing a fast saturable absorber mechanism in fiber lasers via multimodeinterference. We obtain an exact solution for the nonlinear evolution of first and second ordermoments of a laser beam carried by a graded-index multimode fiber, predicting that the spatialself-imaging period does not vary with power. Whereas the amplitude of the oscillation of thebeam width is power-dependent. We have experimentally studied the longitudinal evolution of beamself-imaging by means of femtosecond laser pulse propagation in both the anomalous and the normaldispersion regime of a standard telecom graded-index multimode optical fiber. Light scatteringout of the fiber core via visible fluorescence emission and harmonic wave generation permits us todirectly confirm that the self-imaging period is invariant with power. Spatial shift and splitting ofthe self-imaging process under the action of self-focusing are also emphasized. a r X i v : . [ phy s i c s . op ti c s ] M a y INTRODUCTION
Nonlinear multimode optical fibers (MMFs) are an emerging research field, as they permit new ways for the controlof spatial, temporal and spectral properties of ultrashort pulses of light [1]. As a result, nonlinear MMFs are of interestfor 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 alsoprovide a simply accessible testbed for the study of complex physical phenomena [2]. Among MMFs, graded-indexfibers are of particular interest for nonlinear optics studies, since the reduced modal dispersion leads to relatively longinteraction 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 GRINfibers. SSI leads to periodic oscillations along the fiber of the beam width and intensity. In combination with theKerr effect, SSI leads to a longitudinal periodic modulation of the refractive index of the fiber core, akin to a dynamiclong-period grating. As a result, as recently discussed in a review paper by Agrawal [5], SSI leads to many recentlydiscovered nonlinear effects in GRIN MMFs, such as dispersive wave series emission from multimode femtosecondsolitons [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 themode-locking of fiber lasers, by exploiting multimode interference (MMI) resulting from SSI in a short piece of GRINMMF 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 ( β n − β ) z s = πm n , (1)where β n and β 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 itslength is exactly equal to an integer multiple of z s . However, unless the shape of the input beam exactly matchesthe fundamental mode of the fiber, the input beam is spread in several propagating modes. As a result, the periodicreconstitution 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 evidenceby Zhu et al. [11] by tuning the wavelength of the input signal, and measuring the transmission spectrum of the MMIstructure.Nazemosadat and Mafi proposed to use intensity-dependent differential phase shift among transverse modes, toobtain nonlinear MMI in a short length of GRIN MMF [12]. According to their description, in a nonlinear MMF theself-imaging period should change to z l where ( β n ( I ) − β ( I )) z I = πm n . (2)In a realistic situation, an input laser beam coupled to a GRIN MMF excites several modes, each of them carrying adifferent amount of intensity along the fiber. Therefore the simple description of Eq. (2) should be substituted by thesolution of nonlinear coupled mode equations including self- and cross-phase modulation, as well as four-wave mixingterms [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 observedthat nonlinear MMI acts as a fast saturable absorber (SA) mechanism which permits to obtain mode-locking in highpulse-energy fiber lasers. Specifically, low power signals are strongly attenuated when propagating through the SMSstructure, while high intensity pulses experience a relatively high transmission. Nonlinear MMI-based SAs permitto operate at much higher pulse energies and peak powers than other SA mechanisms, thanks to their high damagethreshold, low cost, simple structure, and mechanical robustness. Nonlinear MMI in a GRIN MMF also permitsspectral 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 demonstratedin recent years [14–22]. In order to test the nonlinear transmission of the SMS structure, the MMF was stretched, andthe SMS transmission was measured as a function of the stretching length, for different input intensities [18]. It wasfound 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 ( (cid:39) ) with respectto the low intensity case ( (cid:39) ). This suggests that nonlinear mode coupling and the resulting beam reshaping couldbe a mode-locking mechanism, rather than the nonlinear variation of the self-imaging period. Indeed, experimentsshow that the excitation of high-order modes (HOMs) in the GRIN fiber is a necessary condition to obtain saturableabsorber action [15, 16]. A study of the impact of the GRIN core diameter (ranging from 20 to 62.5 µm ) revealed thatlarger diameters (hence larger proportions of excited HOMs) increase the output laser power. Another mode-lockingmechanism could result from the longitudinal translation of the self-imaging process, because of self-focusing effects inthe GRIN MMF.Already back in 1992, Karlsson et al. have theoretically studied, by means of the variational method, the dynamicsof 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 thebeam width and its longitudinal phase delay. In this work, we derive an exact solution for the nonlinear evolution offirst and second order moments of a laser beam carried by a GRIN multimode fiber. Our theoretical analysis does notrequire, 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 currentunderstanding of nonlinear MMI based SAs, that the SSI period does not vary with power. Whereas the amplitude ofthe beam width oscillations is power-dependent, which is still sufficient to enable a power-dependent transmission forthe 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 showthat it is possible to work around this problem, by experimentally recording the local nonlinear parametric conversionobtained in a non-collinear geometry, i.e., by using Cherenkov phase matching [24]. As we shall see, the periodicemission of a second harmonic wave, obtained at the core-cladding interfaces, accompanied by multi-photon absorptionleading to wideband fluorescence in the blue spectral domain, provides a clear evidence of the periodic evolution of thepower density in the MMF. We could highlight these processes, we believe for the first time, by using femtosecond laserpulse excitation both in the anomalous and in the normal dispersion regime. This permits us to obtain sufficientlyhigh peak powers, in the multi-MW range, at the points of minimum beam waist of the SSI process. THEORYModel
We consider beam propagation in a multimode optical fiber with a parabolic index profile. The beam dynamics isassumed to be described by a 2D-NLSE with an instantaneous Kerr nonlinearity of the form ∂A∂z − i k (cid:18) ∂ A∂x + ∂ A∂y (cid:19) + i k ∆ ρ ( x + y ) A = i k n n co | A | A, (3)where k = ω n co /c is the wavenumber, ρ is the core radius, ∆ = ( n co − n cl ) / n 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 Ramanscattering and assume that the index profile has an infinite extent (i.e. it does not truncate when the cladding index isreached). Moments of the beam
Now, we show that it is possible to obtain an exact, analytical solution for the moments of the field, whose evolutionis described by the 2D-NLSE (3). In particular, we derive a closed system of equations for the various moments. Themoments may conveniently be obtained by considering them as expectation values of observables and using operatormethods from quantum mechanics. With this machinery the evolution equations for the various moments can be foundin a few lines of calculation using commutator algebra. We introduce a linear Hamiltonian operator ˆ H = 12 k (cid:0) ˆ p x + ˆ p y (cid:1) + β (cid:0) ˆ x + ˆ y (cid:1) − γ ˆ I (4)where β = k ∆ /ρ , γ = k n /n co , ˆ x and ˆ y are transverse position operators, and the nonlinear term ˆ I = | A ( x, y ) | istreated as a potential term that should be determined self-consistently. Here we have also made use of the momentumoperators ˆ p x = − i∂/∂x and ˆ p y = − i∂/∂y . With this pseudo-Hamiltonian we can write Eq. (3) as i ∂A∂z = ˆ HA. (5)Following the procedure described in the Appendix A, we find the sinusoidal solutions for the first-order moments (cid:104) x (cid:105) = a x cos (cid:32)(cid:114) βk z + b x (cid:33) , (cid:104) y (cid:105) = a y cos (cid:32)(cid:114) βk z + b y (cid:33) , (6)where the oscillation period depends on the index parameters, and a x,y , b x,y are integration constants that dependon the initial conditions. Moreover, the radial moment (cid:104) x + y (cid:105) = (cid:104) x (cid:105) + (cid:104) y (cid:105) , associated with the rms-width for acentered beam, satisfies the closed equation d (cid:104) x + y (cid:105) dz = 2 k (cid:104) H (cid:105) − βk (cid:104) x + y (cid:105) . (7)where the Hamiltonian invariant (cid:104) H (cid:105) is defined in Eq. (27). This equation has a solution when (cid:104) H (cid:105) > that describessinusoidal oscillations around a constant average value, viz. (cid:104) x + y (cid:105) ( z ) = 14 β (cid:104) H (cid:105) + a r cos (cid:32)(cid:114) βk z + b r (cid:33) , (8)where a r and b r are integration constants that determine the amplitude and phase of the oscillations. These can befound from the initial beam profile entering the fiber. Whenever the input beam is not in the process of changing,i.e. [ d (cid:104) x + y (cid:105) /dz ] = 0 , one obtains a r = (cid:104) x + y (cid:105) − (cid:104) H (cid:105) / β and b r = 0 . The oscillation amplitude is consequentlya function of the input beam width, while the SSI period is determined by the fiber parameters only, and it remainsequal to its linear value z s = πρ/ √ . It is also seen that the oscillation period of the first-order moments is twicethat of the rms-width.The first integral of Eq. (7) represents a conservation law, M − γk (cid:104) x + y (cid:105)(cid:104) I (cid:105) = const., for the beam qualityfactor M = (cid:104) x + y (cid:105)(cid:104) p x + p y (cid:105) − ( (cid:104) xp x + p x x (cid:105) + (cid:104) yp y + p y y (cid:105) ) , cf. [25]. Furthermore, assuming b r = 0 , we mayrewrite Eq. (8) by using the beam compression parameter (C-parameter) that was introduced by Karlsson et al. [23] (cid:104) x + y (cid:105) ( z ) = (cid:104) x + y (cid:105) (cid:34) (1 + C ) + (1 − C ) cos (cid:32)(cid:114) βk z (cid:33)(cid:35) , C = 12 β (cid:104) H (cid:105)(cid:104) x + y (cid:105) − (9)This parameter measures the relative importance of diffractive broadening and nonlinear self-focusing. The radialmoment is seen to be stationary for C = 1 , while it becomes negative for C < which corresponds to the condition forcollapse. The critical collapse power is reduced in a GRIN fiber with respect to a homogeneous bulk material, becauseof the guiding refractive index profile. In the latter case the Hamiltonian must be negative for collapse to occur for abeam 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 thecorresponding evolution of the associated moment. The normalized simulation parameters are k = 1 . , ∆ = 1 . , ρ =0 . , n co = 2 . and n = 1 . , and the initial field is given by the function A = A sech ( x/a x ) sech ( y/a y ) sin( x ) e i tan − ( y/x ) with coefficients a x = 0 . , a y = 1 . and A = 0 . that has been shifted to the position ( x , y ) = (1 . , . . Thesimulated moment (red dots) is found to be in excellent agreement with the analytical expressions that is shown as ablue 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 Hamiltonianinvariant that is determined by the initial conditions. Assuming a Gaussian initial beam profile of the form A ( x, y,
0) = A exp (cid:34) − (cid:18) xx (cid:19) − (cid:18) yy (cid:19) (cid:35) , (10) FIG. 1. Top panel: Contour slices showing simulated propagation of a complex beam with a vortex charge. The initial beamprofile has been displaced from the center of the fiber (dashed line) and displays periodic oscillations in both position andrms-width; Bottom panel: Evolution of the moment corresponding to the non-centered rms-width the power and Hamiltonian invariants evaluate to (cid:104) P (cid:105) = π | A | x y , (cid:104) H (cid:105) = (cid:104) P (cid:105) (cid:20) k (cid:18) x + 1 y (cid:19) + 2 β ( x + y ) − γ (cid:104) P (cid:105) πx y (cid:21) . (11)The ( /e ) beam width is found by calculating the radial moment (cid:112) (cid:104) x + y (cid:105) / (cid:104) P (cid:105) = (cid:112) ( x + y ) / , which reduces to x = y for a symmetric beam. It is seen that the magnitude of the Hamiltonian is reduced by the nonlinear termfor a focusing nonlinearity. This implies that the average value of the radial beam width will shrink as the power isincreased. In Fig. 2 we illustrate the power-induced variation of the average beam width, i.e. w = (cid:112) (cid:104) H (cid:105) / β (cid:104) P (cid:105) , fordifferent values of the input width x . As can be seen, the average beam width decreases as the power (cid:104) P (cid:105) growslarger: the narrower the input beam size x , 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 beamdiameter of µm (FWHM) so that x = y = 24 . µm , one obtains that associated power-induced variation of theaverage width for the radial moment w occurs on a GW scale (and w shrinks down to zero for (cid:104) P (cid:105) = 2 . GW). Thismeans that the nonlinear reduction of the Hamiltonian is not the root cause of spatial beam cleaning. Note that, inthis case, the average radial beam width approaches the value w = 17 . µm (c.f. core radius ρ = 26 µm ) in the low x = 1.92 mx = 3.84 m x = 5.76 mx = 7.68 m x = 9.60 mx = 11.52 m FIG. 2. Variation of average beam width with power, for unchirped Gaussian initial conditions. power limit. This can be compared with the radial width of the fundamental linear eigenmode solution A ( x, y, z ) = (cid:114) απ exp (cid:20) − α ( x + y ) − i k αz (cid:21) , α = (cid:112) k β, (12)which has w f = (cid:112) (cid:104) x + y (cid:105) / (cid:104) P (cid:105) = 1 / √ α ≈ . µm or . µ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 itreaches zero when (cid:104) P (cid:105) ≈ . MW. However, the beam will not be stationary in the GRIN profile, and oscillationsof the beam width will cause it to collapse at a lower critical power. The moment theory is however exact under theassumption 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. x = 1.92 mx = 3.84 m x = 5.76 mx = 7.68 m x = 9.60 mx = 11.52 m FIG. 3. Variation of beam oscillation amplitude with power, for unchirped Gaussian initial conditions.
To better highlight the nonlinear dependence of the solution for the second-order radial moment Eq. (8), we plot inFig. 3 the input beam power (cid:104) P (cid:105) dependence of the beam width oscillation amplitude A r = 4 βa r / (cid:104) H (cid:105) , for differentvalues of the the input radial width w . As can be seen, if the initial width w = w f , the oscillation amplitude is equalto zero in the linear limit, and it increases with power until it reaches unity for (cid:104) P (cid:105) (cid:39) . M W . This means that thebeam width shrinks to zero, which corresponds to a beam collapse condition. In fact, the C-parameter C = 12 βk x y (cid:18) − γk (cid:104) P (cid:105) x y π ( x + y ) (cid:19) (13)shows that, for a symmetric Gaussian with flat initial phase, the critical collapse power (cid:104) P c (cid:105) = 2 π/ ( γk ) = λ / (2 πn co n ) is independent of the beam width. For a beam with a nonvanishing initial phase-front curvature, our approach alsopermit to obtain the condition for critical collapse to occur. It is interesting to point out that the collapse powerdepends on the incidence angle of the beam into the GRIN MMF. Details about the analytical description of SSI beamevolution 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 beamwidth is larger (smaller) than the width of the fundamental mode w f . A positive (negative) value of A r means that thebeam width is decreasing (increasing) along the fiber with respect to the input value, until it experiences a maximumbeam 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 initialbeam 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 inputpower, we carried out an experimental study of the dynamics of SSI in the nonlinear regime of pulse propagation ina GRIN MMF. By using femtosecond input pulses with peak powers up to the MW power range, we could directlyvisualize the high intensity points inside the MMF, thanks to the associated side scattering of visible higher-harmonicsand fluorescence light. Beamdirection 1 mm 500 µm
FIG. 4. Side scattering imaging of periodic self-imaging in 5 cm of GRIN MMF. From the top, we show the digital microscopepicture, fluorescence at ∼ Anomalous dispersion
We performed two sets of measurements, with two different laser sources. First, we used an ultra-short femtosecondlaser system, involving a hybrid optical parametric amplifier (OPA) of white-light continuum, pumped by a femtosecondYb-based laser, generating 120 fs pulses at 1550 nm, with 25 kHz repetition rate. The input laser beam was focused bya 30 mm focal lens, corresponding to an input beam diameter ( /e ) of 18 µ m, into a 5 cm long multimode standard50/125 GRIN fiber, with relative index difference ∆ = 0 . . FIG. 5. Experimental mesurements of periodic SSI in 5 cm of GRIN MMF, for 5 different input peak power values. Top panels:digital microscope pictures. Bottom panels: relative 1D plot of the fluorescence intensity.
Figures 4 and 5 reveal the presence of multiple peaks of scattered light, corresponding to the nodes of maximumbeam compression in the course of SSI, and highest intensities. This permits us to directly monitor, as shown inFigure 5, the dependence of SSI as a function of the input peak power of the injected pulses. Here the input peakpower 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 showsthat blue fuorescence is periodically scattered from the side of the MMF. The bottom two panels in Figure 4 showthat 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 predictedvalue (i.e., z s = πρ/ √ ≈ µm ). The input pulses generate high-order multimode solitons in the fiber. Thesesolitons undergo, after a propagation distance of about 10 cm, fission into several fundamental solitons under theaction of Raman soliton self-frequency shift and higher-order dispersion. Here we limit ourselves to consider the regimeof multisoliton propagation over the first few cm of GRIN MMF, that is before that soliton fission takes place. Thedynamics 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 ofvisible fluorescence lead to significant nonlinear losses, which induce intensity clamping at the output of the GRINMMF.The invariance of the SSI period with power is well illustrated by Figure 6(a), that shows the input power dependenceof the spatial frequency intensity spectrum of the scattered fluorescence. These results confirm that the SSI periodremains 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, correspondingto 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 ofSection that the ansatz for the transverse mode profile is maintained along the nonlinear beam propagation. On theother 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 axialdirection increases (decreases) along the fiber. The longitudinal dependence of the beam size in the axial directioncould be related with the presence of the significant nonlinear losses, which occur over the first few centimeters of thefiber. To put into evidence the nonlinear transmission properties of the GRIN MMF subject to input peak powers justbelow the critical value for collapse, we reported in Figure 7 the values of the fraction of input coupled energy thatemerges from different lengths of the fiber, as a function of the input peak power. As can be seen, from 1 cm of fiberthe transmission drops below 50% at the highest peak powers close to 5 MW. Whereas for lengths above 7 cm, thetransmission drops below 20% for powers approaching 3 MW. Note that the different curves for fiber lengths above 7cm tend to overlap, indicating that most of the nonlinear loss occurs over the first few centimeters of the fiber. FIG. 6. Experimental measurements of power dependence of (a) self-imaging period; (b) beam transverse and (c) axial dimension:vertical lines represent the variance of the measurements for each input power.FIG. 7. Experimental measurement of nonlinear transmission, for different lengths of the GRIN fiber. Dots: experimental points;curves are provided as a guide to the eye. Normal dispersion
FIG. 8. Experimental results on the non-collinear frequency conversion obtained in a 50/125 GRIN fiber by using a 250 fs lasersource at 1030 nm; (a) schematic representation of the coupling conditions and side image of the fiber, showing periodic emissionof visible light; (b) Spectrum of the emitted light, (c) far field image of the output beam (fiber length: 1 cm, peak power: 2MW).
We also carried out a series of measurements in the normal dispersion regime of the fiber, by using a fiber lasersource at 1030 nm, generating 250 fs pulses with a 30 kHz repetition rate. Here we injected the pulses in a relativelyshort, 1 cm section of a 50/125 µ m GRIN MMF. We obtained at the nodes (that is, at the points of minimum beamwaist 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 arescattered outside the fiber cladding (see Figure 8). The SHG is due to the presence of a weak quadratic nonlinearitybecause of the Ge doping ions.Beyond the observation of the self-imaging periodicity by means of non-collinear frequency conversions, we investigatedpossible longitudinal distortions introduced by a self-focusing regime (see Figure 9). By increasing the input beampower (up to 5.2 MW), and keeping the on-axis excitation, we initiated a self-focusing propagation regime, startingfrom the first node of the self-imaging process.This spatial trapping, which resembles a Townes soliton [27], propagates over hundreds of micrometers beforerecovering its diffracting nature under the effect of nonlinear losses introduced by frequency conversions. Because of itsintensity-dependent nature, this extreme event can drastically modulate the initial mode beating, by introducing botha pulse breaking process and a shift of the self-imaging periodicity. As a result, the intense beam of Figure 9(b) istransformed 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 presenceof SSI, which leads to sharp intensity peaks along the periodic evolution of the field in the GRIN MMF. The measuredperiodicity of the light intensity in the fiber matches well the value obtained with experiments at 1550 nm, as well asthe 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, asit 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 coresection, but is either localized on the central part or on a side. Thus, the longitudinal modulation imprinted on the1 FIG. 9. Experimental side images of a 50/125 optical fiber excited by a 250 fs laser source at 1030 nm showing the periodicself-imaging process; (a) image recorded at 2.1 MW; (b) 5.2 MW; (c) output near and far fields for an input peak power of 2.1MW.FIG. 10. Experimental side images of a 50/125 optical fiber excited by a 250 fs laser source at 1030 nm for 2.1 MW of peakpower; (a) axial coupling; (b) off-axis excitation. propagating light inside the fiber is directly dependent of the initial coupling conditions, i.e. of the combination ofinitially excited modes which can self-modulate because of the mode beating and the Kerr effect, and initiate a spatialexchange of energy between them. This mechanism is at the origin of the spatial self-cleaning effect, which has beenlargely reported in the literature [9, 10]. Additionally, this experiment demonstrates that the transverse position of thenodes of the periodic beating can be transversely tuned by means of the initial coupling conditions. This brings anadditional degree of freedom to realize a nonlinear saturable absorption mechanism, eventually favouring the emergenceof a high-order mode, in order to benefit of the higher dispersion regime when light is coupled back in the multimodefiber. 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 exactsolution for the first and second order moments of a laser beam, describing both the period and the amplitude of2the beam width oscillations along the fiber. The theory also permits to analytically predict the critical power forbeam critical self-focusing, or collapse. We experimentally studied the longitudinal evolution of beam self-imaging bymeans of femtosecond laser pulse propagation in both the anomalous and the normal dispersion regime of a standardGRIN fiber. By observing light scattering out of the fiber core via visible fluorescence emission and harmonic wavegeneration, 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 toall-optical beam processing with multimode fibers.
FUNDING
European Research Council (ERC) (740355); Agence Nationale de la Recherche (ANR) (ANR-18-CE080016-01);CILAS Company (ArianeGroup) by means of the shared X-LAS laboratory; Swedish Research Council (Grant No.2017-05309); Russian Ministry of Science and Education (14.Y26.31.0017);
ACKNOWLEDGMENTS
We thank Fabrizio Frezza for helpful discussions.
DISCLOSURES
The authors declare no conflicts of interest.
Appendix
Momentum operators must satisfy certain commutation relations [ˆ x i , ˆ x j ] = 0 , [ˆ p i , ˆ p j ] = 0 , [ˆ x i , ˆ p j ] = iδ ij , [ˆ x i , ˆ I ( x, y )] = 0 , (14)where i, j denotes combinations of the x, y operators, and the commutator is defined by [ ˆ f , ˆ g ] = ˆ f ˆ g − ˆ g ˆ f . We assumethat (cid:126) = 1 and that the commutators obey the standard Lie algebra [ ˆ f , ˆ g ] = − [ˆ g, ˆ f ] , [ ˆ f , α ˆ g + β ˆ h ] = α [ ˆ f , ˆ g ] + β [ ˆ f , ˆ h ] , [ ˆ f , [ˆ g, ˆ h ]] + [ˆ g, [ˆ h, ˆ f ]] + [ˆ h, [ ˆ f , ˆ g ]] = 0 , (15)together with the relation [ ˆ f , ˆ g ˆ h ] = ˆ g [ ˆ f , ˆ h ] + [ ˆ f , ˆ g ]ˆ h .To derive the moment equations we make use of the fact that the integrals can be expressed as expectation values ofoperators (cid:104) q (cid:105) = (cid:90) A ˆ qA ∗ ds, (16)where ds = dxdy . As a result, the evolution equation for each moment are obtained from d (cid:104) q (cid:105) dz = i (cid:104) [ˆ q, ˆ H ] (cid:105) . (17)In particular we have that (cid:104) x (cid:105) = (cid:90) x | A | ds, (cid:104) p x (cid:105) = i (cid:90) (cid:18) ∂A∂x A ∗ − A ∂A ∗ ∂x (cid:19) ds, (cid:104) I (cid:105) = (cid:90) | A | ds, (18) (cid:104) x (cid:105) = (cid:90) x | A | ds, (cid:104) xp x + p x x (cid:105) = i (cid:90) x (cid:18) ∂A∂x A ∗ − A ∂A ∗ ∂x (cid:19) ds, (cid:104) p x (cid:105) = (cid:90) (cid:12)(cid:12)(cid:12)(cid:12) ∂A∂x (cid:12)(cid:12)(cid:12)(cid:12) ds, (19)and analogous for y .3For the first-order moments we find that d (cid:104) x (cid:105) dz = − k (cid:104) p x (cid:105) , d (cid:104) p x (cid:105) dz = 2 β (cid:104) x (cid:105) , (20) d (cid:104) y (cid:105) dz = − k (cid:104) p y (cid:105) , d (cid:104) p y (cid:105) dz = 2 β (cid:104) y (cid:105) . (21)These equations have sinusoidal solutions (see Eqs. (6))For the second-order moments we have that d (cid:104) x (cid:105) dz = − k (cid:104) xp x + p x x (cid:105) , d (cid:104) y (cid:105) dz = − k (cid:104) yp y + p y y (cid:105) , (22) d (cid:104) p x (cid:105) dz = 2 β (cid:104) xp x + p x x (cid:105) − iγ (cid:104) [ p x , I ] (cid:105) , d (cid:104) p y (cid:105) dz = 2 β (cid:104) yp y + p y y (cid:105) − iγ (cid:104) [ p y , I ] (cid:105) , (23) d (cid:104) xp x + p x x (cid:105) dz = − k (cid:104) p x (cid:105) + 4 β (cid:104) x (cid:105) + γ (cid:104) I (cid:105) , d (cid:104) yp y + p y y (cid:105) dz = − k (cid:104) p y (cid:105) + 4 β (cid:104) y (cid:105) + γ (cid:104) I (cid:105) , (24)where we have used that (cid:104) [ xp x + p x x, I ] (cid:105) = (cid:104) [ yp y + p y y, I ] (cid:105) = i (cid:104) I (cid:105) and it remains to evaluate the commutators betweenthe square momenta and ˆ I . From the definition we have that (cid:104) [ p x , I ] (cid:105) = (cid:90) ∂ ( | A | ) ∂x (cid:18) ∂A∂x A ∗ − A ∂A ∗ ∂x (cid:19) ds, (25)and similar for y . The commutators are related through (here some care is necessary due to the nonlinearity) d (cid:104) I (cid:105) dz = − i k (cid:0) (cid:104) [ p x , I ] (cid:105) + (cid:104) [ p y , I ] (cid:105) (cid:1) . (26)It can be verified that the proper Hamiltonian invariant (cid:104) H (cid:105) = 1 k (cid:0) (cid:104) p x (cid:105) + (cid:104) p y (cid:105) (cid:1) + 2 β (cid:0) (cid:104) x (cid:105) + (cid:104) y (cid:105) (cid:1) − γ (cid:104) I (cid:105) (27)is conserved. In addition to the Hamiltonian, we also have conservation of power and angular momentum, viz. (cid:104) P (cid:105) = (cid:90) | A | ds, (cid:104) yp x − xp y (cid:105) = i (cid:90) (cid:20) y (cid:18) ∂A∂x A ∗ − A ∂A ∗ ∂x (cid:19) − x (cid:18) ∂A∂y A ∗ − A ∂A ∗ ∂y (cid:19)(cid:21) ds. (28) Collapse conditions
The integration constants of Eq. (8) can in general be determined from a r = (cid:18) (cid:104) x + y (cid:105) − β (cid:104) H (cid:105) (cid:19) + k β (cid:18) d (cid:104) x + y (cid:105) dz (cid:19) , (29) tan b r = − (cid:115) k β d (cid:104) x + y (cid:105) dz (cid:30) (cid:18) (cid:104) x + y (cid:105) − β (cid:104) H (cid:105) (cid:19) , (30)with the change in rms-width for the initial beam profile given by d (cid:104) x + y (cid:105) dz = − k [ (cid:104) xp x + p x x (cid:105) + (cid:104) yp y + p y y (cid:105) ] . (31)If the beam has a nonvanishing initial phase-front curvature, b r (cid:54) = 0 , we may rewrite Eq. (9) as (cid:104) x + y (cid:105) ( z ) = (cid:104) x + y (cid:105) b r (cid:34) (1 + C cos b r ) + (1 − C ) cos (cid:32)(cid:114) βk z + b r (cid:33)(cid:35) , (32)4where the C-parameter satisfies C = 14 β (cid:104) H (cid:105)(cid:104) x + y (cid:105) − (cid:104) x + y (cid:105) (cid:115)(cid:18) (cid:104) x + y (cid:105) − β (cid:104) H (cid:105) (cid:19) + k β (cid:18) d (cid:104) x + y (cid:105) dz (cid:19) , (33)and is defined by the condition that the minimum rms-width is (cid:104) x + y (cid:105) = C (cid:104) x + y (cid:105) . Collapse will therefore occurat some point along the beam trajectory when C ≤ . 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 exp (cid:34) − (cid:18) xx (cid:19) − (cid:18) yy (cid:19) (cid:35) . (34)The power and Hamiltonian invariants for this case are obtain as (cid:104) P (cid:105) = Γ (1 / | A | x y , (cid:104) H (cid:105) = 4 π (cid:104) P (cid:105)√ (1 / (cid:20) k (cid:18) x + 1 y (cid:19) + β ( x + y ) − γ (cid:104) P (cid:105) πx y (cid:21) , (35)while the beam width and C-parameter becomes w f = 2 / Γ(1 / (cid:113) π ( x + y ) , C = 32 βk x y (cid:18) − γk (cid:104) P (cid:105) x y π ( x + y ) (cid:19) . (36)Also here it is seen that the critical collapse power is independent of the beam width for a symmetric beam. [1] K. Krupa, A. Tonello, A. Barthélémy, T. Mansuryan, V. Couderc, G. Millot, P. Grelu, D. Modotto, S. A. Babin, andS. Wabnitz, APL Photonics , 110901 (2019), https://doi.org/10.1063/1.5119434.[2] A. Picozzi, G. Millot, and S. Wabnitz, Nat. Photonics , 289 (2015).[3] A. Mafi, Journal of Lightwave Technology , 2803 (2012).[4] H. Talbot, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science , 401 (1836).[5] G. P. Agrawal, Optical Fiber Technology , 309 (2019).[6] L. G. Wright, D. N. Christodoulides, and F. W. Wise, Nat. Photonics , 306 (2015).[7] L. G. Wright, S. Wabnitz, D. N. Christodoulides, and F. W. Wise, Phys. Rev. Lett. , 223902 (2015).[8] K. Krupa, A. Tonello, A. Barthélémy, V. Couderc, B. M. Shalaby, A. Bendahmane, G. Millot, and S. Wabnitz, Phys. Rev.Lett. , 183901 (2016).[9] K. Krupa, A. Tonello, B. M. Shalaby, M. Fabert, A. Barthélémy, G. Millot, S. Wabnitz, and V. Couderc, Nat. Photonics , 234 (2017).[10] L. G. Wright, Z. Liu, D. A. Nolan, M.-J. Li, D. N. Christodoulides, and F. W. Wise, Nat. Photonics , 771 (2016).[11] X. Zhu, A. Schülzgen, H. Li, L. Li, L. Han, J. V. Moloney, and N. Peyghambarian, Opt. Express , 16632 (2008).[12] E. Nazemosadat and A. Mafi, J. Opt. Soc. Am. B , 1357 (2013).[13] F. Poletti and P. Horak, J. Opt. Soc. Am. B , 1645 (2008).[14] S. Fu, G. Shi, Q. Sheng, W. Shi, X. Zhu, J. Yao, R. A. Norwood, and N. Peyghambarian, Opt. Express , 11282 (2016).[15] Z. Wang, D. N. Wang, F. Yang, L. Li, C. Zhao, B. Xu, S. Jin, S. Cao, and Z. Fang, Journal of Lightwave Technology ,5280 (2017).[16] H. Li, Z. Wang, C. Li, J. Zhang, and S. Xu, Opt. Express , 26546 (2017).[17] F. Yang, D. N. Wang, Z. Wang, L. Li, C.-L. Zhao, B. Xu, S. Jin, S.-Y. Cao, and Z.-J. Fang, Opt. Express , 927 (2018).[18] Z. Wang, D. N. Wang, F. Yang, L. Li, C.-L. Zhao, B. Xu, S. Jin, S.-Y. Cao, and Z.-J. Fang, Opt. Lett. , 2078 (2018).[19] U. Teğin and B. Ortaç, Opt. Lett. , 1611 (2018).[20] F. Zhao, Y. Wang, H. Wang, X. Hu, W. Zhang, T. Zhang, and Y. Cai, Laser Physics , 085104 (2018).[21] F. Zhao, Y. Wang, H. Wang, Z. Yan, X. Hu, W. Zhang, T. Zhang, and K. Zhou, Scientific Reports. , 16369 (2018).[22] G. Chen, W. Li, G. Wang, W. Zhang, C. Zeng, and W. Zhao, Photon. Res. , 187 (2019).[23] M. Karlsson, D. Anderson, and M. Desaix, Opt. Lett. , 22 (1992).[24] P. Cerenkov, Dokl, Akad, Nauk, SSSR , 451 (1934).[25] D. Dragoman, Appl. Opt. , 4142 (1996).[26] T. Hansson, D. Anderson, M. Desaix, and M. Lisak, Optics Communications , 3422 (2011).[27] R. Y. Chiao, E. Garmire, and C. H. Townes, Phys. Rev. Lett.13