Optical solitons and vortices in fractional media: A mini-review of recent results

The article produces a brief review of some recent results which predict stable propagation of solitons and solitary vortices in models based on the nonlinear Schroedinger equation including fractional one- or two-dimensional diffraction and cubic or cubic-quintic nonlinear terms, as well as linear potentials. The fractional diffraction is represented by fractional-order spatial derivatives of the Riesz type, defined in terms of the direct and inverse Fourier transform. In this form, it can be realized by spatial-domain light propagation in optical setups with a specially devised combination of mirrors, lenses, and phase masks. The results presented in the article were chiefly obtained in a numerical form. Some analytical findings are included too -- in particular, for fast moving solitons, and results produced by the variational approximation. Also briefly considered are dissipative solitons which are governed by the fractional complex Ginzburg-Landau equation.


I. INTRODUCTION AND THE BASIC MODELS
Nonlinear Schrödinger equations (NLSEs) give rise to soliton families in a great number of realizations [1]- [8], many of which originate in optics. The introduction of fractional calculus in NLSEs has drawn much interest since it was proposedoriginally, in the linear form -as the quantum-mechanical model, derived from the respective Feynman-integral formulation, for particles moving by Lévy flights [9,10]. Experimental implementation of fractional linear Schrödinger equations has been reported in condensed matter [11,12] and photonics [13], in the form of transverse dynamics in optical cavities. Realization of the propagation dynamics of light beams governed by this equation was proposed in Ref. [14], followed by its extension for models including complex potentials subject to the condition of the parity-time (PT ) symmetry [15,16]. Generally, optical settings modeled by linear and nonlinear fractional Schrödinger equations may be considered as a specific form of artificial photonic media.
The objective of this article is to present a short review of models based on fractional NLSEs and some states produced by them. The review is not drafted to be a comprehensive one, the bibliography not being comprehensive either; rather, it selects several recent results which seem interesting in terms of the general soliton dynamics.
The basic form of the one-dimensional (1D) fractional NLSE model, written for amplitude Ψ of the optical field, is Here z and x are the scaled propagation distance and transverse coordinate, g > 0 (g < 0) is the coefficient representing the self-focusing (defocusing) cubic (Kerr) nonlinearity, and V (x) is a trapping potential which may be included in the model. The fractional-diffraction operator with the Lévy index α (in the original fractional Schrödinger equation, it characterizes the hopping motion of the quantum particle [9]) is defined, in the form of the Riesz derivative [45], by means of the juxtaposition of the direct and inverse Fourier transform [9]- [13], [19]: In particular, the above-mentioned concatenation of mirrors, lenses, and phase masks in optical cavities makes it possible to implement the effective fractional diffraction approximated by operator (2) [13]. The Lévy index in Eq. (1) may take values arXiv:2108.11461v1 [nlin.PS] 25 Aug 2021 where α = 2 corresponds to the normal paraxial diffraction, represented by operator −∂ 2 /∂x 2 . Equation (1) with g > 0 (selffocusing) gives rise to the critical collapse at α = 1 (and to the supercritical collapse at α < 1), which destabilizes all possible soliton solutions [23,39], therefore only values α > 1 are usually considered (the linear fractional Schrödinger equation with α = 1 admits analytical solutions based on Airy functions [13]). In the case of quadratic self-focusing nonlinearity, the critical collapse occurs at α = 1/2, hence the Lévy index may take values 1/2 < α ≤ 2 in that case. The quadratic self-focusing term appears on the right-hand side of Eq. (1), in the form of −ε|ψ|ψ with ε > 0, as the correction induced by quantum fluctuations [46] to the 1D fractional Gross-Pitaevskii equation for the gas of Lévy-hopping particles [47].
Besides the real waveguiding potential, V (x), Eq. (1) may include a PT -symmetric complex potential, with even real and odd imaginary parts: where * stands for the complex conjugate. As usual, in terms of optical waveguides the odd imaginary part of potential (4) represents spatially separated and mutually balanced gain and loss elements. The two-dimensional (2D) version of the underlying equation (1), with two transverse coordinates, x and y, is written as In this case, the fractional-diffraction operator is defined as cf. Eq. (2). The rest of the article is organized as follows. In Section II, basic numerical and analytical results are presented for soliton families generated by the 1D fractional NLSE (1) in the free space (without the trapping potential, V = 0). Other essential numerical findings for 1D solitons are presented in Section III. In particular, these are results produced by trapping potentials, and dissipative solitons obtained in the framework of the fractional complex Ginzburg-Landau equation (CGLE) with the cubicquintic nonlinearity. Selected results for 2D solitons, including ones with embedded vorticity and necklace-shaped soliton clusters, are collected in Section IV. To secure stability of the solitons, the 2D model includes the quintic self-defocusing term, in addition to the cubic self-focusing one; the stabilization of 2D solitons can also be provided by the parabolic trapping potential, see Eqs. (47) and (55) below. The paper is concluded by Section V.

II. BASIC SOLITON FAMILIES IN THE ONE-DIMENSIONAL FRACTIONAL MEDIUM
A. Stationary states and analysis of their stability Steady-state solutions to Eq. (1) with propagation constant −µ are looked for in the form of with real function U (x) satisfying the stationary equation, (here, the external potential is dropped, V = 0). The stationary states are characterized by their power (alias norm), Stability of these states was investigated by taking the perturbed solution as where a(x) and b * (x) are components of small complex perturbations, and λ is an instability growth rate (which may be complex). Substituting this ansatz in Eq. (1), one derives the linearized equations for a and b: The underlying stationary solution (7) is stable provided that all eigenvalues λ produced by Eq. (11) have zero real parts.
B. The quasi-local approximation for modes carried by a rapidly oscillating continuous wave To better understand the purport of the fractional-diffraction operator defined by Eq. (2), it is instructive to consider its action on function Ψ(x) built as a slowly varying envelope ψ(x) multiplying a rapidly oscillating continuous-wave carrier, with large wavenumber P . In particular, this ansatz may be used to construct solutions for rapidly moving modes, see Eq. (15) below (in fact, in the spatial domain the "moving" modes are ones tilted in the (x, z) plane). The substitution of wave form (12) in the nonlocal expression on the right-hand side of Eq. (2) readily leads to the following result, in the form of a quasi-local expression expanded in powers of small parameter 1/P : The substitution of expansion (13), truncated at some n = n max , in Eq. (1) and cancellation of the common factor e iP x leads to the familiar NLSE with higher-order dispersion (diffraction) terms, which correspond to n ≥ 3 in Eq. (13). In particular, setting ψ = const demonstrates that the dispersion relation produced by the linearized version of Eq. (8) is in agreement with the analysis of Laskin [10]. Equation (14) implies that the semi-infinite bandgap, with µ < 0, may be populated by solitons, in the framework of the nonlinear equation. As concerns the term ∼ i∂ψ/∂x, corresponding to n = 1, which appears after the substitution of expansion (13) in Eq. (1), it suggests one to rewrite the NLSE in the reference frame moving with the group velocity i.e., to replace x byx in the resulting equation (in fact, the group velocity is the tilt of spatial solitons). Note that, for values of the Lévy index belonging to interval (3), the group velocity is large for large P , while the effective second-order diffraction coefficient, as determined by the term with n = 2 in Eq. (13), is small, The respectively transformed NLSE (1) (without the external potential, V (x) = 0), truncated at n max = 2, is whereψ (x, z) = exp ((i/2)|P | α z) · ψ (x, z). The obvious soliton solution of Eq. (18) with the self-focusing nonlinearity, g > 0, written in terms of norm (9), isψ It follows from Eq. (17) that this soliton is a narrow one, with the width estimated as W ∼ N −1 |P | −(2−α) . Further, in this case the terms representing higher-order diffraction terms with n ≥ 3, which originate from the expansion (13), are relatively small perturbations decaying ∼ |P | −n(α−1) .
C. The scaling relation and variational approximation (VA) for soliton families Equation (8) gives rise to an exact scaling relation between the soliton's power and propagation constant: with a constant N 0 (α) (its particular value is N 0 (α = 2) = 2 √ 2, see Eq. (31) below). The fact that relation (20) satisfies the commonly known Vakhitov-Kolokolov (VK) criterion, dN/dµ < 0 (21) [49,50] at α > 1 suggests that the respective soliton family may be stable, see further details in Fig. 3 below. The case of α = 1, which corresponds to the degenerate form of relation (20), with N (µ) ≡ const, implies the occurrence of the above-mentioned critical collapse, which makes all solitons unstable (cf. the commonly known cubic NLSE in the 2D space with the normal diffraction, α = 2, in which the family of Townes solitons, destabilized by the critical collapse, has a single value of the norm [50]). In the case of α < 1, the solitons generated by Eq. (8) are made strongly unstable by the presence of the supercritical collapse, similar to solitons of the usual cubic NLSE in three dimensions [50]. Localized solutions of the fractional NLSE can be looked for in an approximate analytical form by means of the variational approximation (VA) [23,39,51]. To introduce it, note that Eq. (8) for real U (x), with the fractional diffraction operator defined as per Eq. (2), can be derived from the Lagrangian, The simplest form of the variational ansatz approximating the solitons sought for is based on the Gaussian (cf. Ref. [53]), with real amplitude A, width W , and the power calculated as per Eq. (9), (the calligraphic font denotes quantities pertaining to VA). The substitution of the ansatz in Lagrangian (22) yields the corresponding effective Lagrangian, which may be conveniently written with squared amplitude A 2 replaced by the norm, according to Eq. (24) [39]: where Γ is the Gamma-function. Then, values of N and W are predicted by the Euler-Lagrange equations, Particular examples of the solitons shapes predicted by the VA, and their comparison to the numerically found counterparts are shown in Fig. 1. At a fixed value of the Lévy index, soliton families are characterized by dependences N (µ). An example of such a VA-predicted dependence and its numerical counterpart are displayed in Fig. 2 for α = 1.5, which demonstrates a sufficiently good accuracy of the VA.
An essential result of the VA is the prediction of the fixed norm for the quasi-Townes solitons, for which the critical collapse takes place in the free space at α = 1 [39]: This result is similar to the well-known VA prediction for the norm of the Townes solitons in the 2D NLSE with the cubic self-focusing [53], Reprinted with permission from Ref. [39]. Copyright 2020 Elsevier.
FIG. 2. Dependence of the soliton's norm N on the propagation constant, k ≡ −µ, for the Lévy index α = 1.5 and self-focusing coefficient g = 1, as predicted by the VA based on the Gaussian ansatz (23) and Euler-Lagrange equations (26). The corresponding dependence obtained from the numerical solution of Eq. (8) is included too. Reprinted with permission from Ref. [39]. Copyright 2020 Elsevier. [53], while the respective numerical value is [50], the relative error of the VA being ≈ 7%. It is shown below in Fig. 3(c) that the numerically found counterpart of the variational value (27) is Thus, the relative error of the VA prediction (27) is ≈ 13%. A relatively large size of the error is a consequence of the complicated structure of the equation with the fractional diffraction, especially for values of α which are taken far from the normal-diffraction limit, α = 2. The VA also makes it possible to predict if dependence N (α) for a fixed value of µ is growing or decaying. To make this conclusion, one should take into account that, at α = 2, the usual NLSE with the cubic nonlinearity gives rise to the following N (µ) dependences for the commonly known soliton solutions: Comparing it with the variational value (27), one concludes that the N (α) dependence is growing at −µ > −µ crit ≈ 1/4, and decaying at −µ < −µ crit . This prediction agrees with numerical results displayed below in Fig. 3(c).  (20). The value marked by the arrow is the µ-independent one (27), predicted by the VA for the degenerate family of the quasi-Townes solitons. Its numerically found counterpart is given by Eq. (30). The stability and evolution of the solitons labeled by B1-B4 are displayed in Fig. 4. The plots are borrowed from Ref. [47] (unpublished).

D. Numerical findings
Generic numerically obtained results for soliton families produced by Eqs. (1) and (8) Fig. 3(a). Note that this value of α is taken far from the normal-diffraction limit, α = 2, and close to the collapse boundary, α = 1, in order to demonstrate the setting in which the fractional character of the diffraction is essential. It is seen that the numerically computed dependence exactly follows the analytical prediction given by Eq. (20), with a properly fitted constant N 0 (α). Blue and red segments in the N (µ) curve identify, respectively, stable and unstable soliton subfamilies.
Further, a curve representing a typical dependence N (α) for a fixed propagation constant, −µ = 1.5, which is also split in stable and unstable segments, is exhibited in Fig. 3(c). Unlike the N (µ) dependence, this one cannot be predicted in an exact analytical form. Nevertheless, as mentioned above, the VA based on ansatz (23) predicts the approximate value given by Eq. (27) for the degenerate (µ-independent) norm of the quasi-Townes solitons at α = 1, which is close enough to its numerical counterpart (30). At α = 2, the numerical value N numer (µ = 1.5, α = 2) ≈ 3.46 is in complete agreement with the exact analytical value given by Eq. While the possibility of the collapse in fractional NLSE (1) gives rise to instability of all solitons at α ≤ 1 in the free space (with V = 0), the parabolic trapping potential, helps to create stable solitons even in this case [34]. This finding is similar to the well-known fact that the parabolic trap lifts the norm degeneracy of the Townes solitons and stabilizes their entire family against the critical collapse in the framework of the usual two-dimensional NLSE with the cubic self-attraction [54,55]. Moreover, the same potential makes it possible to predict the existence of stable higher-order (multi-peak) solutions of Eq. (1). Such solutions may be considered as a nonlinear extension of various excited bound states maintained by the parabolic trapping potential in the linear Schrödinger equation. The latter finding is similar to the ability of the 2D parabolic potential to partly stabilize a family of trapped vortex solitons with winding number n = 1 (i.e., the lowest-order excited states in 2D), in the framework of the cubic self-attractive NLSE [54,55]. The higher-order solitons of orders n = 1, 2, 3, ... , supported by Eq. (8) with potential (32), may be approximated by the commonly known stationary wave functions of eigenstates of the quantum-mechanical harmonic oscillator (corresponding to α = 2), where H n (x) are Hermite polynomials, and A is an amplitude. An example of a stable dipole-mode soliton, corresponding to n = 1, and its fit provided by Eq. (33) with properly chosen A, is displayed in Fig. 5 for α = 1, which corresponds to the limit of the collapse-induced instability in the free space. Another noteworthy effect induced by the potential term in the fractional NLSE was recently considered in Ref. [42], which addressed Eq. (1) with a double-well potential, viz., with V 0 > 0 and width w. Equation (1) with this potential may support solitons which are symmetric or antisymmetric with respect to the two potential wells. An effect previously studied in detail in the framework of the usual NLSE (with α = 2) is that, when the soliton's power exceeds a critical value, i.e., the nonlinearity is strong enough, symmetric solitons become unstable and are replaced by stable asymmetric ones, in the case of the self-focusing nonlinearity, i.e., g > 0 in Eq. (1), while antisymmetric solitons remain stable with the increase of their norm [56]. Alternatively, an antisymmetry-breaking bifurcation takes place with antisymmetric solitons (while the symmetric ones remain stable) when the soliton's norm exceeds a critical value in the case of the self-defocusing nonlinearity, corresponding to g < 0 in Eq. (1). Such symmetry-and antisymmetry-breaking bifurcations (phase transitions) in the fractional NLSE with potential (34) are shown, severally, in Figs. 6 and 7 for the Lévy index α = 1.1, which makes the system essentially different from the usual NLSE, corresponding to α = 2. As shown in Fig. 8(a), simulations of the perturbed evolution, performed for the present model in Ref. [42], demonstrate that an unstable symmetric soliton (in the case of the self-attractive nonlinearity) breaks its symmetry and spontaneously transforms into a dynamical state which is close to either one of stable asymmetric solitons existing with the same power (a chiral pair of such solitons is shown in Fig. 6(b)). On the other hand, simulations of the evolution of unstable antisymmetric solitons (in the case of the self-repulsion) produce a different result, as shown in Fig. 8(b): the soliton does not tend to spontaneously transform into either one of the mutually chiral asymmetric solitons, but instead oscillates between them.
Furthermore, the action of the trapping potential (34), even if it is essentially weaker than the parabolic potential (32), provides an effect similar to the above-mentioned one induced by potential (32), viz., stabilization of symmetric, antisymmetric, and asymmetric solitons in the case of α ≤ 1, when all solitons are unstable in the free space, due of the occurrence of the critical (α = 1) or supercritical (α < 1) collapse. In particular, it was demonstrated that the symmetry-breaking bifurcation, quite similar to the one displayed in Fig. 6 (thus, including stable soliton branches) takes place in the same system at α = 0.8 [42].

B. Dissipative solitons produced by the fractional complex Ginzburg-Landau equation (CGLE)
An essential extension of the fractional NLSE was proposed in Ref. [39], in the form of the equation with complex coefficients (including the coefficient in front of the fractional-diffraction term), i.e., the fractional CGLE. It models the propagation of light FIG. 7. (a) Blue, red, and black lines show, respectively, the norm (power, here denoted P ) of antisymmetric, asymmetric, and symmetric solitons vs. their propagations constants (denoted β ≡ −µ here) in the fractional system based on Eq. (1) with the double-well potential (34) and g = −1 (the self-defocusing nonlinearity). Other parameters are the same as in Fig. 6. The solid and dotted segments of the blue line represent, respectively, stable and unstable subfamilies of antisymmetric solitons, below and above the symmetry-breaking bifurcation. The branches of antisymmetric and asymmetric solitons are completely stable. Panels (b), (c), and (d) display, severally, a pair of stable asymmetric solitons (mirror images of each other), an unstable antisymmetric soliton, and a stable symmetric one, all taken at β = 1.12 (these solitons correspond to dots b1,2, c, and d in panel (a)). Reprinted with permission from Ref. [42]. Copyright 2020 Elsevier. in waveguides which, in addition to the fractional diffraction, include losses and gain: where δ > 0, µ > 0, and β ≥ 0 account for, respectively, the linear, quintic, and diffraction losses (fractional diffusion), ε > 0 is the cubic gain, while ν ≥ 0 represents quintic self-defocusing, if it is present in the medium, and the cubic self-focusing term is the same as in Eq. (1) with g = 1.
Equation (35) gives rise to dissipative solitons. Unlike those in the conservative NLSE, dissipative solitons exist not in families parametrized by the propagation constant, but as isolated localized solutions, whose stability is the major issue [39]. Systematic simulations of Eq. (35) have produced charts of different dynamical states, which are displayed in Fig. 9. The charts demonstrate a broad parameter region in which stable dissipative solitons appear.
It is worthy to note that interaction between stable in-phase dissipative solitons, initially separated by some distance, leads to their merger into a single one, as shown in Fig. 10. It is seen that the merger is extremely slow in the case of the usual diffraction (α = 2), while the fractional diffraction essentially accelerates the process, due to the fact that the respective operator (2) actually makes the interaction between the separated solitons nonlocal, i.e., stronger.

A. Stationary vortex solitons
As concerns the propagation of light beams with a 2D transverse structure in fractional media, the transmission of Airy waves, including ones carrying intrinsic vorticity (the optical angular momentum), was analyzed in the framework of the linear variant of Eq. (5) with g = V = 0 [57,58]. In particular, it is worthy to note that the tightest self-focusing of the Airy-ring input is attained at the value of the Lévy index α ≈ 1.4, which is essentially different from α = 2 corresponding to the normal 2D diffraction.
Dynamics of vortex solitons and vorticity-carrying ring-shaped soliton clusters was recently addressed in Refs. [31,35], and [48] in the framework of the fractional NLSE with the cubic-quintic nonlinearity: Unlike the NLSE with the usual diffraction (α = 2), the substitution of ansatz (37) in Eq. (36) with the fractional diffraction operator is not trivial. Nevertheless, using the definition (6) of the operator, the respective calculation can be performed, demonstrating that the vortex ansatz is compatible with Eq. (36) (in other words, the angular-momentum operator commutes with the system's Hamiltonian), the resulting equation being where the fractional radial Laplacian is obtained in a cumbersome form: with Bessel function J 0 . Localized solutions of Eq. (38) are characterized by the 2D norm, In the absence of the self-defocusing quintic term, a straightforward corollary of Eq. (38) with the cubic-only nonlinearity is a scaling relation between the norm and propagation constant, cf. Eq. (20): Obviously, at all values α ≤ 2 (i.e., at all relevant values of the Lévy index), this relation does not satisfy the VK criterion (21), that is why the inclusion of the quintic self-defocusing term is necessary for the stability of the resulting solitons, both fundamental (s = 0) and vortical ones. Furthermore, due to the stabilizing effect of the quintic term, Eq. (36) supports stable localized states even in the case of α ≤ 1, when the cubic self-focusing alone gives rise to the collapse in the 1D setting. This stabilizing effect, which is maintained by the quintic nonlinearity in the free space, may be compared to the above-mentioned stabilization mechanism, provided by the trapping potential. The distribution of the local power in soliton solutions with embedded vorticity (winding number) s features a ring-like shape (see examples below in the left column of Fig. 13). Families of such solutions with s = 1, 2, and 3, produced by numerical FIG. 11. The propagation constant, β ≡ −µ, vs. the norm (power) of the 2D solitons with vorticities s = 1, 2, and 3. The power is defined here as P ≡ 2 2/α N2D, see Eq. (40). The results are produced by the numerical solution of Eq. (38) for fixed values of the Lévy index: α = 1 (a), 1.3 (b), 1.5 (c), 1.7 (d), 1.9 (e), and 2 (f). The latter case, corresponding to the usual two-dimensional NLSE, is included for comparison with the results produced by the fractional diffraction. Note that values of the power are shown on the logarithmic scale, which is different in different panels. Blue and red segments represent stable and unstable vortex states, respectively. Reprinted with permission from Ref. [31]. Copyright 2020 Elsevier. solution of Eq. (38), are represented by dependences of the 2D norm on the propagation constant, which are displayed in Fig.  11. These dependences also highlight relatively small stable subfamilies of the vortex solitons. The stability was identified by the computation of eigenvalues, using linearized equations for small perturbations, and verified by direct simulations of the perturbed evolution [31].
In the limit of the diverging norm, the propagation constant in all panels of Fig. 11 attains the commonly known limit value for the NLSE with the cubic-quintic nonlinearity, β max = 3/16 [59], which does not depend on the spatial dimension. Further, it is seen in the figure that the soliton families exist at values of the norm (power) exceeding a certain threshold value, although the stability boundary corresponds to much larger values of the norm. Actually, this feature (i.e., the non-existence of solitons with the norm falling below a finite threshold value) is well known in models where the self-focusing nonlinear term, if acting alone, gives rise to the supercritical collapse, such as the three-dimensional NLSE with the normal diffraction [64]. Indeed, in collapse-free models solitons with a vanishingly small norm have a vanishing amplitude and diverging size, the latter fact implying µ → −0. However, Eq. (41) demonstrates that, at any α < 2, the norm of the solitons is diverging, rather than vanishing, at µ → −0, hence the limit of N → 0 cannot be attained. Dependences of the threshold norm on the Lévy index for s = 1, 2, and 3 are displayed in Fig. 12. In the limit of α = 2, i.e., in the case of the cubic-quintic NLSE with the normal 2D diffraction, the respective threshold values N As well as their fundamental counterparts, families of vortex Townes solitons are degenerate, in the sense that their norm takes the single value, which depends on s but does not depend on the propagation constant. Actually, these values may be approximated by an analytical formula obtained in Ref. [67]: N (s) thr (α = 2) ≈ 4 √ 3πs. In simulations of Eq. (36) those vortex solitons which are unstable split in sets of separating fragments, as shown in Fig. 13. This is the usual instability-development scenario for vortex solitons in NLSEs with diverse nonlinearities [7].  42)), vs. the Lévy index α, for vorticities s = 1, 2, and 3. In the limit of α = 2, which corresponds to the normal (non-fractional) diffraction, the threshold values N (s) thr (α = 2) coincide with the norms of the Townes solitons with the embedded vorticity, see Eq. (43). Reprinted with permission from Ref. [31]. Copyright 2020 Elsevier.

B. Dynamics of vortical clusters
In addition to the axisymmetric stationary vortex solitons, it is relevant to consider the dynamics of necklace-shaped clusters. At z = 0 they are composed as sets (circular chains) of M fundamental solitons (with zero intrinsic winding numbers), with equal distances between them, and global vorticity, S = 1, 2, 3, ..., imprinted onto the cluster: Here U 0 (|r − r m |) is the stationary shape of the fundamental soliton with the center placed at point r = r m , and R is the radius of the cluster. Previously, patterns of the necklace type have drawn much interest in studies of models based on NLSEs with normal diffraction [60][61][62]. In particular, robust soliton necklaces supported by the cubic-quintic nonlinearity (the same as in Eq. (36)) were addressed in Ref. [63].
The formation and stability of the circular soliton chains in the framework of Eq. (36) was studied in detail in Ref. [35]. Parameters of the input, taken in the form of Eqs. (45) and (46), were chosen so that the axisymmetric vortex soliton with the same vorticity S, and the norm equal to the total norm of the initial soliton chain, is stable, in terms of Fig. 11. A set of typical results produced by this input is displayed in Fig. 14. It is seen in row (a) of the figure that, in the absence of the overall vorticity (S = 0), the input set of solitons promptly fuses into a stable fundamental soliton. On the other hand, if the vorticity is too large, S ≥ 2, no bound state is formed, and the initial soliton chain quickly expands, as observed in Fig. 14(b).
The formation of quasi-stable slowly rotating necklaces is possible in the case of S = 1. In Fig. 14(c), the necklace state features slow rotation in combination with periodic oscillations in the radial direction. On the other hand, Fig. 14(d) demonstrates an example of a rotating necklace with a nearly permanent shape, which remains stable in the course of very long propagation. Indeed, z = 1000 in this case is tantamount, roughly, to ∼ 50 diffraction (Rayleigh) lengths of the input pattern, which, in turn, may be estimated as z Rayleigh ∼ (2R) α . C. Stabilization of 2D solitons by the trapping potential As said above, all soliton solutions of Eq. (5) with the cubic self-focusing nonlinearity (g = 1) in the free space (V = 0) are completely unstable at all values of the Lévy index, α ≤ 2. Instead of the quintic self-defocusing term, stabilization may be provided by the 2D parabolic potential, cf. Eq. (32). In an analytical form, this possibility can be demonstrated for small values of Ω 2 in Eq. (47), competing with small fractality parameter (2 − α), which determines the proximity to the normal diffraction (α = 2). First, the effect of the weak trapping potential can be taken into regard by means of the respective VA for the NLSE including this potential, the regular diffraction (α = 2), and cubic self-focusing [68]. The VA for the stationary wave function with zero vorticity (see Eq. (37) with s = 0) is based on the 2D Gaussian ansatz, with the norm (see Eq. (40)) cf. Eqs. (23) and (24). The respective expression for the effective Lagrangian is cf. Eq. (25). The Euler-Lagrange equations following from here, ∂L eff /∂ (W, N ) = 0, lead to an expression showing how potential (47) with small Ω 2 lifts the norm degeneracy of the 2D Townes solitons: Here, the first term corresponds to the above-mentioned VA prediction for the degenerate norm of the Townes solitons in the 2D cubic NLSE with the normal diffraction (α = 2), see Eq. (28). In the lowest approximation, a small effect of Ω 2 may be neglected in the corresponding VA-predicted relation between the soliton's width W and propagation constant −µ : While Eq. (51) pertains to α = 2, the expansion of the free-space (V = 0) scaling relation (41) for small fractality, 0 < 2 − α 1, yields where the value of N for the Townes soliton at α = 2 is substituted by the VA value (28). Combining Eqs. (51) and (53) yields an expression which makes it possible to analyze the competition of the stabilizing and destabilizing effects produced, respectively, by the trapping potential and the weak fractality: The stable part of the family of weakly-fractional solitons characterized by dependence (54) may be identified by means of the VK criterion (21). It predicts that stable solitons are those with the propagation constant −µ > 0 subject the following constraint: Thus, according to Eq. (52), stable solitons should be wide enough, W 2 > √ 2 − α/ (2Ω). This conclusion is natural, as the soliton should be sufficiently wide to feel the stabilizing effect of the trapping potential.

V. CONCLUSION
Studies of the wave propagation in fractional media have made remarkable progress, starting from the fractional linear Schrödinger equation, which was introduced by Laskin [9] for quantum-mechanical particles moving by Lévy flights. The fractality in that equation is represented by the derivative of the Riesz type, characterized by the Lévy index, α. Actually, this is an integral operator (2), defined by means of the combination of direct and inverse Fourier transforms. This operator may be approximately reduced to a combination of usual local derivatives when it acts on the wave function built as a rapidly oscillating continuous-wave carrier multiplied by a slowly varying envelope, as shown by Eq. (13).
The next step was the realization of the fractional Schrödinger equation as one governing the paraxial wave propagation in optical setups emulating the fractional diffraction [13]. The implementation of the fractional Schrödinger equation in optics has made it natural to include the Kerr (as well as non-Kerr) nonlinearities, thus arriving at the fractional NLSE. The nonlinear equations render it possible to predict various self-trapped modes, such as solitons and solitary vortices, in these settings. The present article offers a brief review of some recent theoretical results on this topic, while experimental observations of fractional solitons have not been published, as yet. The modes considered in the review include basic 1D solitons and 2D solitary vortices in the free space, supported by the cubic and cubic-quintic nonlinearities, respectively. The 1D solitons are considered in the interval of values of the Lévy index 1 < α ≤ 2 (the fractional NLSE with the cubic self-focusing term in 1D gives rise to the collapse at α ≤ 1, while α = 2 corresponds to the usual diffraction term). All 2D solitons supported by cubic nonlinearity in the free space are unstable at α ≤ 2.
Also considered are 1D solitons in external potentials, which make it possible to produce stable higher-order (multiple-peak) solitons and study the spontaneous symmetry breaking in double-well potentials. In addition to that, trapping potentials may stabilize 1D solitons at α < 1 and 2D ones at α < 2. For the 2D setting, dynamics of necklace-shaped soliton clusters carrying overall vorticity is included too. Proceeding to the fractional cubic-quintic CGLE, the review briefly surveys 1D dissipative solitons predicted by that equation.
This mini-review is far from being a comprehensive one. Among other topics related to the fractional NLSEs are, inter alia, PT -symmetric solitons maintained by complex potentials subject to condition (4). In particular, breaking and restoration of the solitons' PT symmetry was recently considered in Ref. [16]. An interesting possibility is to introduce fractional equations of the NLSE type with quadratic or quadratic-cubic nonlinearity [47,69]. Also theoretically considered were various modes supported by a modulationally stable background field (continuous wave), such as dark solitons and vortices, "bubbles", and W-shaped solitons [51,70].
The theory may be developed in other directions. In particular, NLSEs of a different type, that account for temporal propagation of optical waves under the action of the fractional group-velocity dispersion, were developed in Refs. [71] and [72]. The analysis has produced solitons in those models as well.