Spatiotemporal dissipative solitons and vortices in a multi-transverse-mode fiber laser

We introduce a model for spatiotemporal modelocking in multimode fiber lasers, which is based on the (3+1)-dimensional cubic-quintic complex Ginzburg-Landau equation (cGLE) with conservative and dissipative nonlinearities and a 2-dimensional transverse trapping potential. Systematic numerical analysis reveals a variety of stable nonlinear modes, including stable fundamental solitons and breathers, as well as solitary vortices with winding number $n=1$, while vortices with $n=2$ are unstable, splitting into persistently rotating bound states of two unitary vortices. A characteristic feature of the system is bistability between the fundamental and vortex spatiotemporal solitons.

Mode-locking effects and formation of dissipative solitons have been traditionally (and successfully) modeled using various realizations of the complex Ginzburg-Landau equation (cGLE), its ubiquitous version being one with the cubic-quintic nonlinearity [1,28]. The quintic term supports the subcritical transition to the lasing state, thereby providing conditions for mode-locking and the creation of stable solitons at subcritical pumping levels [28]. Locking of the longitudinal modes in single-transverse-mode fiber lasers is described by (1+1)-dimensional cGLE, while the consideration of multi-mode fibers opens the way to the realm of the (3+1)D cGLE. The lateral confinement of the fiber modes immediately requires inclusion of the trapping potential acting in two transverse dimensions, which is a basic parabolic potential for GRIN fibers. This potential is a feature creating a new playground relative to the quite extensive studies of "light bullets" [29][30][31][32] and vortex solitons [33][34][35][36][37] in (3+1)D cGLE. Below we use the (3+1)D cGLE including conservative and dissipative cubic-quintic terms and a 2D parabolic potential.
Measuring the propagation distance along the fiber, beam's width, and time in units of suitably defined scales L, w and τ , the cGLE is written as (1) Here z, t and x, y are the dimensionless distance, time, and transverse coordinate, and time is defined in the reference frame moving with the pulse's group velocity. Further, k is the GVD coefficient, β the spectral-filtering coefficient, k 0 the propagation constant of the fundamental mode, g the GRIN parameter, and ε ≡ Lε > 0 the dimensionless linear loss over length L. If ψ is the dimensionless field envelope, then the scaled (dimensionless) parameters η = Lη , α = Lα > 0, ν = Lν and µ = Lµ > 0 determine ratios between L and the nonlinear lengths associated with the Kerr effect, cubic gain, quintic nonlinearity, and nonlinear saturation, respectively. We now fix L = 1/(k 0 gw 2 ), w 2 = L/k 0 and τ 2 = Lk to write the spatiotemporal cGLE (1) with the cubic-quintic terms in the scaled form: Nontrivial solutions can be supported by the cubic gain in the competition with the linear and quintic losses, if the gain coefficient exceeds a threshold value [39] α = Aα thr , α thr ≡ (2 √ εµ) with A > 1. We here fix values A = 1.5 and β = ν = 0, which adequately represent the generic case, while the linear loss and nonlinear gain, ε and µ, are varied.
Thus, in this work we focus on predicting STML scenarios in the setup based on the interplay of the cubic-quintic dissipative nonlinearity with the anomalous GVD, transverse diffraction, parabolic GRIN potential, and Kerr self-focusing, in the absence of temporal and spatial filtering. In fact, the prediction of patterns which are stable without the help of the filtering is an essential novelty of the analysis.

II. NUMERICAL RESULTS
We look for stationary localized solutions, in the form of ψ(x, y, t, z) = φ(x, y, t) exp (iqz), where q is a real propagation constant, and function φ(x, y, t) satisfies the equation following from the substitution of this ansatz in Eq. (2) with β = ν = 0: Splitting complex equation (4) into real and imaginary parts, solutions for localized states have been produced by means of the modified squared-operator method [40]. The solutions are characterized by the value of their total energy, i.e., integral norm,

A. Fundamental dissipative solitons
To determine numerical solutions and the propagation constant simultaneously, the Gaussian input was used, φ 0 = 1.5 sech (x 2 + y 2 ) exp(−1/2t 2 ). The existence and stability of numerically found fundamental (zero-vorticity) dissipative solitons, produced by Eq. (4), is summarized in Fig. 1(a), which displays the numerically found q and P (integral power) vs. the linear-loss parameter, ε, for a fixed value of the quintic-loss coefficient, µ = 1, while the cubic gain is fixed as per Eq.  The breathers are robust dynamical modes, which keep integrity in the course of indefinitely long subsequent evolution. Their stability regions in the parameter plane of (µ, ε) are also drawn in Fig. 1(b). Note, in particular, that the region of destruction of all modes (i.e., onset of spatial "turbulence" as a result of the simulated evolution) occupies, quite naturally, an area in which the loss and gain coefficients [recall that the latter one is taken as per Eq.
(3)] are too large to maintain the stability of regular modes, such as solitons and breathers.
The same peculiarity is observed in stability charts displayed below for vortex states in Figs. Localized vortex states with winding number (topological charge) n = 1 were constructed as solutions of Eq. (4), starting with input, written in polar coordinates (r, θ) (in the plane of (x, y)) as φ 0 = 2.5r sech (r 2 ) sech (t) exp (iθ). As a result, two coexisting vortex families, one unstable and one stable, have been found. Figure 5(a) displays the respective numerically found values of q and integral power P versus ε. In particular, the stable (red) branch was found in the interval of 0 < ε < 0.4 for µ = 1.0. Unstable vortex states split into fundamental solitons (in the blue segment), or transform into chaotically oscillating modes (magenta) or moving breathers, which keep the initial vorticity (yellow), or suffer destruction (black).

5(b) and 10(b).
The results are summarized in Fig. 5(b) in terms of the (µ, ε) plane, where the shaded area The stability chart for vortices with n = 1 in the plane of (µ, ε). The central shaded area in (b) is populated by chaotically oscillating localized modes, see an example in Fig. 8. Note that the stability area of vortices is completely covered by the region in which the fundamental solitons are stable [cf. Fig. 1(b)], i.e., it is a bistability area.

III. CONCLUSIONS
In this work, we have proposed and investigated a model for the passively mode-locked laser, based on the graded-index nonlinear multimode fiber. The model is based on the (3 + 1)D cGLE (complex Ginzburg-Landau equation) with the complex cubic-quintic nonlinearity and the transverse two-dimensional trapping potential, which represents the graded-index structure. A novelty in comparison with previously studied models is that stable stationary devices is expected to assist the development of methods to optimize and tune temporal, spatial and spectral characteristics of short spatiotemporal pulses generated in multimode fiber lasers. [1] P. Grelu, N. Akhmediev, "Dissipative solitons for mode-locked lasers," Nature Photonics 6(2),