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; Stefan Wabnitz; Vincent Couderc

doi:10.1364/OE.398531

1. 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 discussed in a review paper by Agrawal [5], SSI leads to many recently discovered nonlinear effects in GRIN MMFs, such as dispersive wave sideband series emission from multimode femtosecond solitons [6,7], geometric parametric instability of continuous wave (CW) beams in the normal dispersion regime [8,9], and spatial beam self-cleaning phenomena [10–12]. In addition, SSI has been recently widely exploited for the mode-locking of fiber lasers, by means of 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

(1)$$(\beta_n-\beta_1)z_s=\pi m_n,$$

where $\beta _n$ and $\beta _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 a superposition of modes with different propagation constants, which affects the output beam divergence. The presence of self-imaging in a MMF was indirectly put into evidence by Zhu et al. [13] 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 [14]. According to their description, in a nonlinear MMF the self-imaging period should change to $z_I$ where

(2)$$(\beta_n(I)-\beta_1(I))z_I=\pi 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 [15]. Equivalently, one may use a two-dimensional nonlinear Schrödinger equation (2D-NLSE) [10].

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 [16–24]. 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 [20]. 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 ($\simeq 75\%$) with respect to the low intensity case ($\simeq 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 [17,18]. A study of the impact of the GRIN core diameter (ranging from 20 to 62.5 $\mu 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 a variational approach (VA), the dynamics of self-imaging in a GRIN MMF under the combined action of diffraction, nonlinearity, and parabolic index profile [25]. They derived an approximate analytical solution for a multimode beam, which permits to describe the evolution of the beam width, and its longitudinal phase delay. The VA is based on the assumption that the beam keeps the initial shape (e.g., Gaussian) along the whole propagation distance. Notable predictions of the VA are that the SSI period $z_s$ is independent of power, whereas the amplitude of the beam intensity and width oscillations is power-dependent. An important consequence of the invariance of $z_s$ is that the position of GPI sidebands remains a constant, as the beam power is increased [8,26]. In this work, we advance the theory of SSI in GRIN MMFs, by deriving an exact solution for the nonlinear evolution of first and second order moments of a laser beam of arbitrary shape. Specifically, our theoretical analysis does not require, unlike methods based on the VA, that the beam maintains a specific shape throughout its propagation in the MMF: it applies even to the case of speckled beams. This permits us to confirm, in full generality, and in contrast with current understanding of fiber laser mode-locking based on nonlinear MMI, 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 from a MMI device, in order to act as equivalent fast SA. Remarkably, our method also permits us to derive the power threshold for catastrophic self-focusing in a GRIN MMF for a beam of arbitrary incidence.

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 photo-luminescence side-scattered by the GRIN MMF. A physical mechanism for light scattering is the emission of a second harmonic (SH) wave, obtained at the core-cladding interfaces, leading to nonlinear parametric conversion in a non-collinear geometry, i.e., by using Cherenkov phase matching [27]. The generation of a SH is accompanied by multi-photon absorption by defects of the Germanium-doped fiber [28], leading to wideband photo-luminescence in the blue spectral domain. The side-scattering of blue light at the points of peak intensity provides a clear evidence of the periodic evolution of the beam propagating 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 intensities, in the multi-MW range, at the points of minimum beam waist of the SSI process.

2. Theory

2.1 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

(3)$$\frac{\partial A}{\partial z} - i\frac{1}{2k_0}\left(\frac{\partial^2 A}{\partial x^2}+\frac{\partial^2 A}{\partial y^2}\right) + i\frac{k_0\Delta}{\rho^2}(x^2+y^2)A = i\frac{k_0n_2}{n_{co}}|A|^2A,$$

where $k_0 = \omega _0n_{co}/c$ is the wavenumber, $\rho$ is the core radius, $\Delta = (n_{co}^2-n_{cl}^2)/2n_{co}^2$ 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).

2.2 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 (see Appendix A for definitions) 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

(4)$$\hat{H} = \frac{1}{2k_0}\left(\hat{p}_x^2+\hat{p}_y^2\right) + \beta\left(\hat{x}^2+\hat{y}^2\right)-\gamma \hat{I}$$

where $\beta = k_0\Delta /\rho ^2$, $\gamma = k_0n_2/n_{co}$, $\hat {x}$ and $\hat {y}$ are transverse position operators, and the $z$-dependent nonlinear operator $\hat {I} = |A(x,y,z)|^2$ is treated as a potential term that should be determined self-consistently. Here we have also made use of the momentum operators $\hat {p}_x = -i\partial /\partial x$ and $\hat {p}_y = -i\partial /\partial y$. With this pseudo-Hamiltonian we can write Eq. (3) as

(5)$$i\frac{\partial A}{\partial z} = \hat{H}A.$$

Following the procedure described in the Appendix, we find the sinusoidal solutions for the first-order moments

(6)$$\langle x\rangle = a_x\cos{\left(\sqrt{\frac{2\beta}{k_0}}z+b_x\right)}, \qquad \langle y\rangle = a_y\cos{\left(\sqrt{\frac{2\beta}{k_0}}z+b_y\right)},$$

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 $\langle x^2+y^2\rangle = \langle x^2\rangle + \langle y^2\rangle$, associated with the root mean square (RMS) width for a centered beam, satisfies the closed equation

(7)$$\frac{d^2\langle x^2+y^2\rangle}{dz^2} = \frac{2}{k_0}\langle H_0\rangle - \frac{8\beta}{k_0}\langle x^2+y^2\rangle.$$

where the Hamiltonian invariant $\langle H_0\rangle$ is defined in Eq. (28) of the Appendix. This equation has a solution when $\langle H_0\rangle > 0$ that describes sinusoidal oscillations around a constant average value, viz.

(8)$$\langle x^2+y^2\rangle (z) = \frac{1}{4\beta}\langle H_0\rangle + a_r\cos\left(\sqrt{\frac{8\beta}{k_0}}z+b_r\right),$$

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\langle x^2+y^2\rangle /dz]_0 = 0$, one obtains $a_r = \langle x^2+y^2\rangle _0 - \langle H_0\rangle /4\beta$ 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=\pi \rho /\sqrt {2\Delta }$. 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 - \gamma k_0\langle x^2+y^2\rangle \langle I\rangle = \textrm {const.}$, for the beam quality factor $M^2 = \langle x^2+y^2\rangle \langle p_x^2+p_y^2\rangle -\frac {1}{4}\left (\langle xp_x+p_xx\rangle +\langle yp_y+p_yy\rangle \right )^2$, cf. [29]. 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. [25]

(9)$$\langle x^2+y^2\rangle (z) = \frac{\langle x^2+y^2\rangle_0}{2}\left[\left(1+C\right)+\left(1-C\right)\cos\left(\sqrt{\frac{8\beta}{k_0}}z\right)\right], \qquad C = \frac{1}{2\beta}\frac{\langle H_0\rangle}{\langle x^2+y^2\rangle_0}-1.$$

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 [30].

2.3 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 radial moment. The normalized simulation parameters are $k_0 = 1.3, \Delta = 1.7, \rho = 0.8, n_{co} = 2.0$ and $n_2 = 1.5$, and the initial field is given by the function $A = A_0~\textrm {sech}(x/a_x)\textrm {sech}(y/a_y)\sin (x)e^{i2\textrm {tan}^{-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).

Fig. 1. Top panel: Contour slices showing simulated propagation of a complex beam with a vortex charge. The initial beam profile has been displaced from the center of the fiber (dashed line) and displays periodic oscillations in both position and width; Bottom panel: Evolution of the moment corresponding to the non-centered RMS-width.

Download Full Size | PDF

2.4 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

(10)$$A(x,y,0) = A_0\exp\left[-\frac{1}{2}\left(\frac{x}{x_0}\right)^2-\frac{1}{2}\left(\frac{y}{y_0}\right)^2\right],$$

the power and Hamiltonian invariants evaluate to

(11)$$\langle P_0\rangle = \pi|A_0|^2x_0y_0, \qquad \langle H_0\rangle = \frac{\langle P_0\rangle}{2}\left[\frac{1}{k_0}\left(\frac{1}{x_0^2}+\frac{1}{y_0^2}\right)+2\beta(x_0^2+y_0^2)-\frac{\gamma \langle P_0\rangle}{\pi x_0y_0}\right].$$

The RMS-width is found by calculating the power normalized radial moment $\sqrt {\langle x^2+y^2\rangle /\langle P_0\rangle } = \sqrt {(x_0^2+y_0^2)/2}$, which reduces to the ($1/e$) beam width $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 = \sqrt {\langle H_0\rangle /4\beta \langle P_0\rangle }$, for different values of the input width $x_0$. As can be seen, the average beam width decreases as the power $\langle P_0\rangle$ grows larger: the narrower the input beam size $x_0$, the faster the nonlinear decrease of the beam width with power.

Fig. 2. Variation of average beam width with power, for unchirped Gaussian initial conditions.

Download Full Size | PDF

When using fiber parameters from the numerical simulations reported in Ref. [10], and assuming a fixed input beam diameter of $x_{\textrm {FWHM}} = 40~\mu m$ full width at half maximum (FWHM) so that $x_0 = y_0 = x_{\textrm {FWHM}}/(2\sqrt {\log {2}}) = 24.0~\mu m$, one obtains that the 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 $\langle P_0\rangle =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~\mu m$ (c.f. core radius $\rho = 26~\mu m$) in the low power limit. This can be compared with the radial width of the fundamental linear eigenmode solution

(12)$$A(x,y,z) = \sqrt{\frac{\alpha}{\pi}}\exp\left[-\frac{1}{2}\alpha(x^2+y^2)-i\frac{1}{k_0}\alpha z\right], \qquad \alpha = \sqrt{2k_0\beta},$$

which has $x_f=\sqrt {\langle x^2+y^2\rangle /\langle P_0\rangle } = 1/\sqrt {\alpha } \approx 5.76~\mu m$ or $9.59~\mu m$ FWHM (see corresponding orange curve in Fig. 2). Note that, for an input beam perfectly matching the fundamental mode shape, the average beam width shrinks until it reaches zero when $\langle P_0\rangle \approx 16.55~\textrm {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 $\langle P_0\rangle$ dependence of the beam width oscillation amplitude $A_r = 4\beta a_r/\langle H_0\rangle$, for different values of the input radial width $x_0$. As can be seen, if the initial width $x_0=x_f$, the oscillation amplitude is equal to zero in the linear limit, and it increases with power until it reaches unity for $\langle P_0\rangle \simeq \mathrm{8.3}\: MW$. This means that the beam width shrinks to zero, which corresponds to a beam collapse condition. In fact, the C-parameter

(13)$$C = \frac{1}{2\beta k_0x_0^2y_0^2}\left(1-\frac{\gamma k_0\langle P_0\rangle x_0y_0}{\pi(x_0^2+y_0^2)}\right)$$

shows that, for a symmetric Gaussian with flat initial phase, the critical collapse power $\langle P_c\rangle = 2\pi /(\gamma k_0)=\lambda ^2/(2\pi n_{co}n_2)$ is independent of the beam width. For a beam with a nonvanishing initial phase-front curvature, our approach also permits 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 are provided in the Appendix.

Fig. 3. Variation of beam oscillation amplitude with power, for unchirped Gaussian initial conditions.

Download Full Size | PDF

Figure 3 also shows that, in the low-power limit, the oscillation amplitude is positive (negative) when the input beam width $x_0$ is larger (smaller) than the width of the fundamental mode $x_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 analytical description to the case of a super-Gaussian initial beam profile: corresponding results are reported in Appendix A.2.

3. 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 photo-luminescence light.

3.1 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 70-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 $\mu$m, into a 5 cm long multimode standard 50/125 GRIN fiber, with relative index difference $\Delta =0.0102$.

Figures 4 and 5 reveal the presence of multiple peaks of scattered light, corresponding to the nodes of maximum beam compression in the course of SSI, and highest intensities. This permits us to directly monitor, as shown in Fig. 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.

Fig. 4. Side scattering imaging of periodic self-imaging in 5 cm of GRIN MMF. From the top, we show the digital microscope picture, photo-luminescence at $\sim$7 MW of input peak power of 120 fs pulses, the side scattering of blue/violet radiation, and its zoom.

Download Full Size | PDF

Fig. 5. Experimental measurements of periodic SSI in 5 cm of GRIN MMF, for five different input peak power values of 120 fs pulses. Top panels: digital microscope pictures. Bottom panels: relative 1D plot of the photo-luminescence intensity.

Download Full Size | PDF

Figure 6(a) shows the multiple peaks of scattered light, as well as the spectra at the fiber output for different input peak powers (see Fig. 6(b)), exhibiting a series of dispersive wave sidebands generated from soliton oscillations [6,7], and the spectrum of side-scattered light (Fig. 6(c)). In all spectra it is possible to note the presence of the emission of a broadband photo-luminescence, with a peak around 360 nm.

Fig. 6. Experimental results on non-collinear frequency conversion obtained in a 50/125 GRIN fiber by using a 70 fs laser source at 1550 nm; (a) side image of the fiber, showing periodic emission of blue light; (b) output spectra for several input peak powers; (c) spectrum of side-scattered light.

Download Full Size | PDF

Figure 4 shows, in the top panel, the picture taken by a digital microscope. The second panel from the top shows that blue photo-luminescence is periodically scattered from the side of the MMF. The bottom two panels in Fig. 4 show that green light is also scattered by the defects of the fiber cladding.

As can be seen in Fig. 5, the SSI period remains unchanged, and remarkably close to the theoretically predicted value (i.e., $z_s=\pi \rho /\sqrt {2\Delta } \approx 550\; \mu 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 photo-luminescence 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 Fig. 7(a), that shows the input power dependence of the spatial period of the scattered photo-luminescence. These results confirm that the SSI period remains a constant (within experimental measurement errors) even at the highest input peak powers.

Fig. 7. Experimental measurements of power dependence of (a) self-imaging period; (b) luminescence transverse diameter and (c) axial dimension at FWHM: vertical lines represent the variance of the measurements for each input power.

Download Full Size | PDF

We also estimated the power-dependence of the beam size at the points of maximum beam compression, corresponding to the bright spots in Fig. 5. Figure 7(b) reports the diameter of the luminescence signal in the transverse section of the fiber both at FWHM and at $1/e$: as can be seen the beam size does not show a significant dependence on input power, keeping a slightly smaller size than the $18~\mu m$ input beam diameter. This confirms the hypothesis of Section 2 that the ansatz for the transverse mode profile is maintained along the nonlinear beam propagation. On the other hand, Fig. 7(c) shows that the beam size at FWHM in the axial dimension exhibits a power-dependence.

Specifically, Fig. 7(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 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 Fig. 8 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.

Fig. 8. Experimental measurement of nonlinear transmission, for different lengths of the GRIN fiber. Dots: experimental points; curves are provided as a guide to the eye.

Download Full Size | PDF

3.2 Normal dispersion

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 $\mu$m GRIN MMF. We obtained at the nodes (that is, at the points of minimum beam waist and maximum intensity) of SSI a sufficient intensity to trigger, in addition to multi-photon luminescence, also non-collinear (i.e., using Cherenkov phase matching [27]) second-harmonic generation (SHG), both of which are scattered outside the fiber cladding (see Fig. 9). The SHG is due to the presence of a weak quadratic nonlinearity because of the Ge doping ions.

Fig. 9. Experimental results on non-collinear frequency conversion obtained in a 50/125 GRIN fiber by using a 250 fs laser source at 1030 nm; (a) schematic representation of the coupling conditions and side image of the fiber, showing periodic emission of visible light; (b) Spectrum of the emitted light, (c) far field image of the output beam (fiber length: 1 cm, peak power: 2 MW).

Download Full Size | PDF

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 Fig. 10). 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.

Fig. 10. Experimental side images of a 50/125 optical fiber excited by a 250 fs laser source at 1030 nm showing the periodic self-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.1 MW.

Download Full Size | PDF

This spatial trapping, which resembles a Townes soliton [31], 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 Fig. 10(b) is transformed in a double peak of intensity, both oscillating with the same initial period of the self-imaging process.

As shown by Fig. 11, 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.

Fig. 11. Experimental side images of a 50/125 optical fiber excited by a 250 fs laser source at 1030 nm for 2.1 MW of peak power; (a) axial coupling; (b) off-axis excitation.

Download Full Size | PDF

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 Fig. 11, 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 [10,12]. 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 when light is coupled back in the multimode fiber. In this sense, temporal mode-locking on a high-order transverse mode should be possible.

4. 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 photo-luminescence 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.

A. Appendix

Momentum operators must satisfy certain commutation relations

(14)$$[\hat{x}_i,\hat{x}_j] = 0, \qquad [\hat{p}_i,\hat{p}_j] = 0, \qquad [\hat{x}_i,\hat{p}_j] = i\delta_{ij}, \qquad [\hat{x}_i,\hat{I}(x,y)] = 0,$$

where $i,j$ denotes combinations of the $x,y$ operators, and the commutator is defined by $[\hat {f},\hat {g}] = \hat {f}\hat {g}-\hat {g}\hat {f}$. We assume that $\hbar = 1$ and that the commutators obey the standard Lie algebra

(15)$$[\hat{f},\hat{g}] ={-}[\hat{g},\hat{f}], \qquad [\hat{f},\alpha \hat{g}+\beta \hat{h}] = \alpha [\hat{f},\hat{g}]+\beta [\hat{f},\hat{h}], \qquad [\hat{f},[\hat{g},\hat{h}]]+[\hat{g},[\hat{h},\hat{f}]]+[\hat{h},[\hat{f},\hat{g}]] = 0,$$

together with the relation $[\hat {f},\hat {g}\hat {h}] = \hat {g}[\hat {f},\hat {h}]+[\hat {f},\hat {g}]\hat {h}$.

To derive the moment equations we make use of the fact that the integrals can be expressed as expectation values of operators

(16)$$\langle q\rangle = \int A\hat{q}A^*~ds,$$

where $ds = dx dy$, and we omit the power normalization for simplicity. As a result, the evolution equation for each moment are obtained from Heisenberg’s equation of motion

(17)$$\frac{d\langle q\rangle}{dz} = i\langle[\hat{q},\hat{H}]\rangle + {\bigg \langle}\frac{\partial q}{\partial z}{\bigg \rangle}. $$

Where the last term accounts for any explicit z dependence of the operator. In particular we have that

(18)$$\langle x\rangle = \int x|A|^2~ds, \qquad \langle p_x\rangle = \frac{i}{2}\int\left(\frac{\partial A}{\partial x}A^*-A\frac{\partial A^*}{\partial x}\right)~ds, \qquad \langle I(z)\rangle = \int |A|^4~ds,$$

(19)$$\langle x^2\rangle = \int x^2|A|^2~ds, \qquad \langle xp_x+p_xx\rangle = i\int x\left(\frac{\partial A}{\partial x}A^*-A\frac{\partial A^*}{\partial x}\right)~ds, \qquad \langle p_x^2\rangle = \int {\bigg |}\frac{\partial A}{\partial x}{\bigg |}^2~ds,$$

and analogous for $y$.

The zeroth-order moment corresponding to the unity operator is a constant of motion that expresses the conservation of power and is denoted by

(20)$$\langle P_0\rangle = \int |A|^2~ds.$$

For the first-order moments we find that

(21)$$\begin{aligned}& \frac{d\langle x\rangle}{dz} ={-}\frac{1}{k_0}\langle p_x\rangle, & \frac{d\langle p_x\rangle}{dz} = 2\beta\langle x\rangle, \end{aligned}$$

(22)$$\begin{aligned}& \frac{d\langle y\rangle}{dz} ={-}\frac{1}{k_0}\langle p_y\rangle, & \frac{d\langle p_y\rangle}{dz} = 2\beta\langle y\rangle. \end{aligned}$$

These equations have sinusoidal solutions (see Eqs. (6)).

For the second-order moments we have that

(23)$$\begin{aligned}& \frac{d\langle x^2\rangle}{dz} ={-}\frac{1}{k_0}\langle xp_x+p_xx\rangle, & \frac{d\langle y^2\rangle}{dz} ={-}\frac{1}{k_0}\langle yp_y+p_yy\rangle, \end{aligned}$$

(24)$$\begin{aligned}& \frac{d\langle p_x^2\rangle}{dz} = 2\beta\langle xp_x+p_xx\rangle - i\gamma\langle [p_x^2,I]\rangle, & \frac{d\langle p_y^2\rangle}{dz} = 2\beta\langle yp_y+p_yy\rangle - i\gamma\langle [p_y^2,I]\rangle, \end{aligned}$$

(25)$$\begin{aligned}& \frac{d\langle xp_x+p_xx\rangle}{dz} ={-}\frac{2}{k_0}\langle p_x^2\rangle + 4\beta\langle x^2\rangle + \gamma\langle I\rangle, & \frac{d\langle yp_y+p_yy\rangle}{dz} ={-}\frac{2}{k_0}\langle p_y^2\rangle + 4\beta\langle y^2\rangle + \gamma\langle I\rangle, \end{aligned}$$

where we have used that $\langle [xp_x+p_xx,I]\rangle = \langle [yp_y+p_yy,I]\rangle = i\langle I\rangle$ and it remains to evaluate the commutators between the square momenta and $\hat {I}$. From the definition we have that

(26)$$\langle [p_x^2,I]\rangle = \int\frac{\partial (|A|^2)}{\partial x}\left(\frac{\partial A}{\partial x}A^*-A\frac{\partial A^*}{\partial x}\right)~ds,$$

and similar for $y$. The commutators are related through

(27)$$\frac{d\langle I\rangle}{dz} ={-}i\frac{1}{k_0}\left(\langle [p_x^2,I]\rangle+\langle [p_y^2,I]\rangle\right).$$

It can be verified that the proper Hamiltonian invariant

(28)$$\langle H_0\rangle = \frac{1}{k_0}\left(\langle p_x^2\rangle+\langle p_y^2\rangle\right) + 2\beta\left(\langle x^2\rangle+\langle y^2\rangle\right)-\gamma\langle I\rangle$$

is conserved. In addition to power and Hamiltonian, we also have conservation of angular momentum, viz.

(29)$$\langle yp_x-xp_y\rangle = \frac{i}{2}\int\left[y\left(\frac{\partial A}{\partial x}A^*-A\frac{\partial A^*}{\partial x}\right)-x\left(\frac{\partial A}{\partial y}A^*-A\frac{\partial A^*}{\partial y}\right)\right]~ds.$$

A.1. Collapse conditions

The integration constants of Eq. (8) can in general be determined from

(30)$$a_r^2 = \left(\langle x^2+y^2\rangle_0 - \frac{1}{4\beta}\langle H_0\rangle\right)^2 + \frac{k_0}{8\beta}\left(\frac{d\langle x^2+y^2\rangle_0}{dz}\right)^2,$$

(31)$$\tan{b_r} ={-}\sqrt{\frac{k_0}{8\beta}}\frac{d\langle x^2+y^2\rangle_0}{dz}{\bigg /}\left(\langle x^2+y^2\rangle_0-\frac{1}{4\beta}\langle H_0\rangle\right),$$

with the change in the radial moment for the initial beam profile given by

(32)$$\frac{d\langle x^2+y^2\rangle_0}{dz} ={-}\frac{1}{k_0}\left[\langle xp_x+p_xx\rangle_0+\langle yp_y+p_yy\rangle_0\right].$$

If the beam has a nonvanishing initial phase-front curvature, $b_r \neq 0$, we may rewrite Eq. (9) as

(33)$$\langle x^2+y^2\rangle (z) = \frac{\langle x^2+y^2\rangle_0}{1+\cos{b_r}}\left[\left(1+C\cos{b_r}\right)+\left(1-C\right)\cos\left(\sqrt{\frac{8\beta}{k_0}}z+b_r\right)\right],$$

where the C-parameter satisfies

(34)$$C = \frac{1}{4\beta}\frac{\langle H_0\rangle}{\langle x^2+y^2\rangle_0} - \frac{1}{\langle x^2+y^2\rangle_0}\sqrt{\left(\langle x^2+y^2\rangle_0 - \frac{1}{4\beta}\langle H_0\rangle\right)^2 + \frac{k_0}{8\beta}\left(\frac{d\langle x^2+y^2\rangle_0}{dz}\right)^2},$$

and is defined by the condition that the minimum beam width is $\langle x^2+y^2\rangle = C\langle x^2+y^2\rangle _0$. Collapse will therefore occur at some point along the beam trajectory when $C \leq 0$.

A.2. Supergaussian beam profile

To investigate the influence of the beam shape on the dynamics we also consider a supergaussian beam profile

(35)$$A(x,y,0) = A_0\exp\left[-\frac{1}{2}\left(\frac{x}{x_0}\right)^4-\frac{1}{2}\left(\frac{y}{y_0}\right)^4\right].$$

The power and Hamiltonian invariants for this case are obtain as

(36)$$\langle P_0\rangle = \frac{\Gamma^2(1/4)}{4}|A_0|^2x_0y_0, \qquad \langle H_0\rangle = \frac{4\pi\langle P_0\rangle}{\sqrt{2}\Gamma^2(1/4)}\left[\frac{3}{2k_0}\left(\frac{1}{x_0^2}+\frac{1}{y_0^2}\right)+\beta(x_0^2+y_0^2)-\frac{\gamma \langle P_0\rangle}{\pi x_0y_0}\right],$$

while the beam RMS-width and C-parameter becomes

(37)$$x_{\textrm{RMS}} = \frac{2^{1/4}}{\Gamma(1/4)}\sqrt{\pi(x_0^2+y_0^2)}, \qquad C = \frac{3}{2\beta k_0x_0^2y_0^2}\left(1-\frac{2\gamma k_0\langle P_0\rangle x_0y_0}{3\pi(x_0^2+y_0^2)}\right).$$

It is seen that the critical collapse power $\langle P_c\rangle = 3\pi /(\gamma k_0)$ is independent of the beam width also for a symmetric supergaussian beam, and is a factor $3/2$ larger than in the case of an ordinary Gaussian beam.

Funding

European Research Council (740355); Agence Nationale de la Recherche (ANR-18-CE080016-01); CILAS Company (ArianeGroup); Vetenskapsrådet (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.

References

1. K. Krupa, A. Tonello, A. Barthélémy, T. Mansuryan, V. Couderc, G. Millot, P. Grelu, D. Modotto, S. A. Babin, and S. Wabnitz, “Multimode nonlinear fiber optics, a spatiotemporal avenue,” APL Photonics 4(11), 110901 (2019). [CrossRef]

2. A. Picozzi, G. Millot, and S. Wabnitz, “Nonlinear virtues of multimode fibre,” Nat. Photonics 9(5), 289–291 (2015). [CrossRef]

3. A. Mafi, “Pulse propagation in a short nonlinear graded-index multimode optical fiber,” J. Lightwave Technol. 30(17), 2803–2811 (2012). [CrossRef]

4. H. Talbot, “LXXVI. Facts relating to optical science. No. IV,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 9(56), 401–407 (1836). [CrossRef]

5. G. P. Agrawal, “Invite paper: Self-imaging in multimode graded-index fibers and its impact on the nonlinear phenomena,” Opt. Fiber Technol. 50, 309–316 (2019). [CrossRef]

6. L. G. Wright, D. N. Christodoulides, and F. W. Wise, “Controllable spatiotemporal nonlinear effects in multimode fibres,” Nat. Photonics 9(5), 306–310 (2015). [CrossRef]

7. L. G. Wright, S. Wabnitz, D. N. Christodoulides, and F. W. Wise, “Ultrabroadband dispersive radiation by spatiotemporal oscillation of multimode waves,” Phys. Rev. Lett. 115(22), 223902 (2015). [CrossRef]

8. S. Longhi, “Modulational instability and space time dynamics in nonlinear parabolic-index optical fibers,” Opt. Lett. 28(23), 2363–2365 (2003). [CrossRef]

9. K. Krupa, A. Tonello, A. Barthélémy, V. Couderc, B. M. Shalaby, A. Bendahmane, G. Millot, and S. Wabnitz, “Observation of geometric parametric instability induced by the periodic spatial self-imaging of multimode waves,” Phys. Rev. Lett. 116(18), 183901 (2016). [CrossRef]

10. K. Krupa, A. Tonello, B. M. Shalaby, M. Fabert, A. Barthélémy, G. Millot, S. Wabnitz, and V. Couderc, “Spatial beam self-cleaning in multimode fibres,” Nat. Photonics 11(4), 237–241 (2017). [CrossRef]

11. Z. Liu, L. G. Wright, D. N. Christodoulides, and F. W. Wise, “Kerr self-cleaning of femtosecond-pulsed beams in graded-index multimode fiber,” Opt. Lett. 41(16), 3675–3678 (2016). [CrossRef]

12. L. G. Wright, Z. Liu, D. A. Nolan, M.-J. Li, D. N. Christodoulides, and F. W. Wise, “Self-organized instability in graded-index multimode fibres,” Nat. Photonics 10(12), 771–776 (2016). [CrossRef]

13. X. Zhu, A. Schülzgen, H. Li, L. Li, L. Han, J. V. Moloney, and N. Peyghambarian, “Detailed investigation of self-imaging in large core multimode optical fibers for application in fiber lasers and amplifiers,” Opt. Express 16(21), 16632–16645 (2008). [CrossRef]

14. E. Nazemosadat and A. Mafi, “Nonlinear multimodal interference and saturable absorption using a short graded-index multimode optical fiber,” J. Opt. Soc. Am. B 30(5), 1357–1367 (2013). [CrossRef]

15. F. Poletti and P. Horak, “Description of ultrashort pulse propagation in multimode optical fibers,” J. Opt. Soc. Am. B 25(10), 1645–1654 (2008). [CrossRef]

16. S. Fu, G. Shi, Q. Sheng, W. Shi, X. Zhu, J. Yao, R. A. Norwood, and N. Peyghambarian, “Dual-wavelength fiber laser operating above 2 μm based on cascaded single-mode-multimode-single-mode fiber structures,” Opt. Express 24(11), 11282–11289 (2016). [CrossRef]

17. Z. Wang, D. N. Wang, F. Yang, L. Li, C. Zhao, B. Xu, S. Jin, S. Cao, and Z. Fang, “Er-doped mode-locked fiber laser with a hybrid structure of a step-index-graded-index multimode fiber as the saturable absorber,” J. Lightwave Technol. 35(24), 5280–5285 (2017). [CrossRef]

18. H. Li, Z. Wang, C. Li, J. Zhang, and S. Xu, “Mode-locked Tm fiber laser using SMF-SIMF-GIMF-SMF fiber structure as a saturable absorber,” Opt. Express 25(22), 26546–26553 (2017). [CrossRef]

19. F. Yang, D. N. Wang, Z. Wang, L. Li, C.-L. Zhao, B. Xu, S. Jin, S.-Y. Cao, and Z.-J. Fang, “Saturable absorber based on a single mode fiber-graded index fiber-single mode fiber structure with inner micro-cavity,” Opt. Express 26(2), 927–934 (2018). [CrossRef]

20. Z. Wang, D. N. Wang, F. Yang, L. Li, C.-L. Zhao, B. Xu, S. Jin, S.-Y. Cao, and Z.-J. Fang, “Stretched graded-index multimode optical fiber as a saturable absorber for erbium-doped fiber laser mode locking,” Opt. Lett. 43(9), 2078–2081 (2018). [CrossRef]

21. U. Teğin and B. Ortaç, “All-fiber all-normal-dispersion femtosecond laser with a nonlinear multimodal interference-based saturable absorber,” Opt. Lett. 43(7), 1611–1614 (2018). [CrossRef]

22. F. Zhao, Y. Wang, H. Wang, X. Hu, W. Zhang, T. Zhang, and Y. Cai, “High-energy solitons generation with a nonlinear multimode interference-based saturable absorber,” Laser Phys. 28(8), 085104 (2018). [CrossRef]

23. F. Zhao, Y. Wang, H. Wang, Z. Yan, X. Hu, W. Zhang, T. Zhang, and K. Zhou, “Ultrafast soliton and stretched-pulse switchable mode-locked fiber laser with hybrid structure of multimode fiber based saturable absorber,” Sci. Rep. 8(1), 16369 (2018). [CrossRef]

24. G. Chen, W. Li, G. Wang, W. Zhang, C. Zeng, and W. Zhao, “Generation of coexisting high-energy pulses in a mode-locked all-fiber laser with a nonlinear multimodal interference technique,” Photonics Res. 7(2), 187–192 (2019). [CrossRef]

25. M. Karlsson, D. Anderson, and M. Desaix, “Dynamics of self-focusing and self-phase modulation in a parabolic index optical fiber,” Opt. Lett. 17(1), 22–24 (1992). [CrossRef]

26. H. E. L. Aviles, F. O. Wu, Z. Eznaveh, M. A. Eftekhar, F. Wise, R. Correa, and D. N. Christodoulides, “A systematic analysis of parametric instabilities in nonlinear parabolic multimode fibers,” APL Photonics 4(2), 022803 (2019). [CrossRef]

27. P. Cerenkov, “Visible glow under exposure of gamma radiation,” Dokl, Akad, Nauk, SSSR 2, 451 (1934).

28. P. G. Kazansky, H. Inouye, T. Mitsuyu, K. Miura, J. Qiu, K. Hirao, and F. Starrost, “Anomalous anisotropic light scattering in Ge-doped silica glass,” Phys. Rev. Lett. 82(10), 2199–2202 (1999). [CrossRef]

29. D. Dragoman, “Wigner distribution function in nonlinear optics,” Appl. Opt. 35(21), 4142–4146 (1996). [CrossRef]

30. T. Hansson, D. Anderson, M. Desaix, and M. Lisak, “Self-similar collapse and blow up of laser beams in nonlinear kerr media: revisited,” Opt. Commun. 284(13), 3422–3427 (2011). [CrossRef]

31. R. Y. Chiao, E. Garmire, and C. H. Townes, “Self-trapping of optical beams,” Phys. Rev. Lett. 13(15), 479–482 (1964). [CrossRef]

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

Abstract

1. Introduction

2. Theory

2.1 Model

2.2 Moments of the beam

2.3 Vortex beam

2.4 Gaussian beam

3. Experiments

3.1 Anomalous dispersion

3.2 Normal dispersion

4. Conclusion

A. Appendix

A.1. Collapse conditions

A.2. Supergaussian beam profile

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (11)

Equations (37)

Optics Express