Basic fractional nonlinear-wave models and solitons

This review article provides a concise summary of one- and two-dimensional models for the propagation of linear and nonlinear waves in fractional media. The basic models, which originate from fractional quantum mechanics and more experimentally relevant setups emulating fractional diffraction in optics, are based on the Riesz definition of fractional derivatives, which are characterized by the respective Levy indices. Basic species of one-dimensional solitons, produced by the fractional models which include cubic or quadratic nonlinear terms, are outlined too. In particular, it is demonstrated that the variational approximation is relevant in many cases. A summary of the recently demonstrated experimental realization of the fractional group-velocity dispersion in fiber lasers is also presented.

The formal calculus based on the concept of fractional derivatives is a branch of formal mathematics that has been known in the course of ca.200 years.Ca. 20 years ago, this concept had appeared in physics as an ingredient of fractional quantum mechanics, which is based on the Schrödinger equation with the fractional operator of kinetic energy, for the wave function of particles which, in their classical form, move by random Lévy flights.The fractionality of the Schrödinger equations of this type is characterized by the Lévy index (LI) α, which takes values 0 < α ≤ 2. It determines the form of the fractional kinetic-energy operator, as −∂ 2 /∂x 2 α/2 , the usual (nonfractional) quantum mechanics corresponding to α = 2.While fractional quantum mechanics remains far from experimental implementation, much interest has been drawn to the more recent proposal to emulate the fractional diffraction, modeled by the same equations, by means of the wave propagation in specially devised optical cavities.Many theoretical results have been reported for solitons, vortices and other modes supported by the respective fractional Schrödinger equations which include cubic or quadratic optical nonlinearities.Very recently, the first experimental implementation of the effective fractional group-velocity dispersion has been reported for the light propagation in a fiber-laser cavity.The great deal of the current interest to these theoretical and experimental studies suggest relevance of providing a summary of the state of art in the area.Such a concise summary is offered by this article, which outlines basic models of the fractional wave propagation, and a brief overview of basic species of solitons produced by these models.

I. INTRODUCTION: THE CONCEPT OF FRACTIONAL-ORDER DERIVATIVES
Ordinary derivatives are the central concept of the classical calculus, which has been the mathematical basis of modern science in the course of the last 350 years.In addition to the commonly known definition of derivatives of integer orders, mathematical curiosity suggested attempts to introduce derivatives of fractional orders.For example, the Darrieus-Landau instability of plane combustion fronts, propagating in a premixed gas, is determined by the linear dispersion relation between the instability growth rate γ and wavenumber k of one-dimensional (1D) perturbations corrugating the flat combustion front, with a constant C > 0 [1].In Eq. (1, |k| implies that the same instability gain is produced by right-and left-traveling perturbations.In the framework of the phenomenological model [2], one can introduce a model equation for a real order parameter u (x, t) which corresponds to the dispersion relation (1), where −∂ 2 /∂x 2 is an operator that gives rise to factor |k| in the Fourier space, where k is the wavenumber conjugate to the spatial coordinate x running along the flat combustion front.Similarly, 2D perturbations with wave vector k = (k x , k y ) give rise to the instability dispersion relation which suggests to introduce the respective 2D phenomenological equation, where ∆ is the usual 2D Laplacian acting on coordinates (x, y) on the flat combustion front.It is shown in detail below that, while there are different formal definitions of fractional differential operators, such as the one √ −∆ in Eq. ( 4), the definitions which are relevant for the realization in physics are defined in terms of the action of the operators in the Fourier space, similar to the relation between Eqs. ( 2) and (1) in 1D, or between Eqs. ( 4) and (3) in 2D.
Much earlier than fractional derivatives have appeared in physical models, they were considered as a formal generalization of the classical calculus.In a systematic form, the concept of the fractional calculus was elaborated by Niels Henrik Abel in 1823 [3] and Joseph Liouville in 1832 [7].The development of these concepts has led to the abstract definition of what is known as the Caputo derivative of non-integer order α [8,9], where n ≡ [α] + 1, [α] stands for the integer part of the fractional order, {α} ≡ α − [α] is its fractional part, Γ is the Gamma-function, and f (n) is the usual derivative.Note that this definition represents the fractional derivative as integral (nonlocal) operator, rather than as a differential one.As shown below, the fractional derivatives which appear in models originating in physics are defined differently, namely, as Riesz derivatives [10], see Eqs. ( 9) and ( 14) in Section 2. This paper aims to produce a concise overview of linear and nonlinear dynamics in physically relevant models of fractional media.Also included is a summary of the first experimental realization of the concept of the fractional derivatives, in the form of a fiber cavity which emulates fractional group-velocity dispersion (GVD) [11].The presentation is not comprehensive; in particular, while the main objective is to introduce models of linear and nonlinear fractional media which are relevant to physics, the summary of soliton solutions amounts to the presentation of a few examples, many other cases being only briefly mentioned.An earlier review of theoretical results for solitons in fractional models is provided in Ref. [12].

II. THE ADVENT OF FRACTIONAL CALCULUS TO PHYSICS: FRACTIONAL QUANTUM MECHANICS
For the first time, fractional derivatives had appeared in the context of physically relevant models as the mathematical basis of fractional quantum mechanics.This theory was introduced by Laskin [13,14] for nonrelativistic particles which move, at the classical level, by Lévy flights (random leaps).

A. Classical particles moving by Lévy flights
The starting point for the development of fractional quantum mechanics is to consider particles whose classical stochastic motion is performed as a chain of random leaps ("Lévy flights"), rather than in the form of the usual Brownian random walk.In this regime, the mean distance |x| of the particle, which moves by 1D "flights" from the initial position, x = 0, grows with time t as where the Lévy index (LI) α takes values 0 < α ≤ 2 (7) [16].The limit value, α = 2, corresponds to the usual Brownian random walk, while at α < 2 Eq. ( 6) demonstrates that the Lévy flights give rise to faster growth of |x| at t → ∞.As concerns realizations of the Lévy flights in nature, a usual reference is to the life style of sharks in the ocean, which search for food moving by random fast swims, even if sharks are not appropriate objects for considering their motion in the quantum regime, which is the next step of the derivation, leading to the fractional Schrödinger equation (FSE) B. The fractional linear Schrödinger equations for Lévy-flying particles and the fractional Riesz derivative The derivation of an effective Schrödinger equation as a result of the quantization of a particle which moves, at the classical level, by Lévy flight was elaborated by Laskin [13,14] (see also Ref. [15]).It is based on the fundamental formalism which defines quantum mechanics by means of the Feynman's path integration.This formalism represents the quantum dynamics of a particle as a result of the superposition of virtual motions along all randomly chosen trajectories (paths).The superposition is defined as the integral in the space of all paths, ∼ exp [iS(path)] d(path), where S is the classical action corresponding to a particular path [17].To apply this concept, Laskin considered the superposition of the paths which correspond not to the Brownian random walks, but to the Lévy flights.In the 1D setting, the result is FSE for wave function Ψ (x, t).In the scaled form, it is written as where V (x) is the same external potential which appears in the usual Schrödinger equation.
A fundamentally important result of the derivation of Eq. ( 8) by means of the path-integral formalism is that the usual kinetic-energy operator, −(1/2)∂ 2 /∂x 2 , which corresponds to α = 2 in Eq. ( 8), is replaced, at α < 2, by the fractional Riesz derivative [10], which is defined as the juxtaposition of the direct and inverse Fourier transforms (x → p → x) of the wave function, This definition implies that the action of the fractional derivative, −∂ 2 /∂x 2 α/2 , in the Fourier space amounts to the straightforward multiplication: Ψ(p) → |p| α Ψ(p), where is the Fourier transform of the wave function.In definition (9), α is the same LI which determines the classical Lévy flights as per Eq. ( 6).Thus, the derivation of the physically relevant model, viz., FSE, leads to the relatively simple Riesz fractional derivative.As concerns "more sophisticated" varieties of fractional derivatives, defined on abstract mathematical grounds, such as the Caputo derivative (5), they do not emerge naturally in these physical contexts.
In the 2D space, the same derivation gives rise to the 2D FSE for wave function Ψ (x, y, t)) [14]: The fractional operator of the 2D kinetic energy is represented, in the 2D Fourier space (p, q), by multiplier p 2 + q 2 α/2 , i.e., its action amounts to where the 2D Fourier transform of the wave function cf. Eq. (10).Therefore, the action of the 2D operator in the coordinate space is defined as the juxtaposition of the direct and inverse Fourier transforms, with the multiplication as per Eq. ( 12) inserted between them: cf. the 1D equation (9).Stationary solutions to Eq. ( 11), with real energy eigenvalue µ (in terms of Bose-Einstein condensates (BECs), it is the chemical potential ) are looked for in the usual form, Ψ (x, y, t) = e −iµt U (x, y) , (15) with real function U satisfying the stationary equation, or its 1D version, without coordinate y.It is relevant to mention that in the case of potential in the form of the harmonic oscillator, V (x, y) = Ω 2 /2 x 2 + y 2 , the Fourier transform (13) satisfies the Fourier transform of Eq. (11), which takes the form of the usual Schrödinger equation, but written in the Fourier space [4]: with the effective potential (1/2) p 2 + q 2 α/2 , or the 1D version of this equation, in the absence of variable q.In particular, a solution of the 1D version of Eq. ( 17) with α = 1, i.e., with the effective Fourier-space potential (1/2)|p|, can be expressed in terms of the Airy function [4].
The fractional calculus has found another important application to physics in the form of fractional kinetics.This concept, which was elaborated by Zaslavsky et al. [5,6], addresses, in particular, the kinetic theory, based on the fractional Fokker-Planck equations, for dynamics which may be intermediate between integrable and purely chaotic.It predicts fundamental effects such as anomalous transport and superdiffusion, and models real kinetic phenomena such as advection of particles in diverse physical settings.However, this topic is not considered in detail in the present article.
C. A conjecture: fractional Gross-Pitaevskii equations (FGPEs) It is commonly known that the dynamics of BECs in ultracold atomic gases is very accurately approximated, in the framework of the mean-field theory, by the Gross-Pitaevskii equation, in the form of the Schrödinger equation for single-particle wave function ψ, supplemented by the cubic terms which represent the effect of inter-particle collisions [41]: Here m is the atomic mass, and a s is the scattering length which determines two-particle collisions (a s > 0 and a s < 0 correspond to repulsive and attractive interactions between the particles, respectively).An intriguing possibility is to consider a gas of bosonic particles obeying the FSE (8) or its 2D version (11).While a systematic derivation of the respective fractional Gross-Pitaevskii equation (FGPE), with the usual kinetic-energy operator in Eq. ( 18) replaced by its fractional counterpart defined as per Eqs.( 8) or (11), remains a challenging objective, a natural expectation is that the equation sought for will take the following form, in the scaled form: where σ = +1 or −1 corresponds to the self-repulsive or attractive condensate, respectively [12].
A recently elaborated example of the FGPE system is a model of a binary (two-component) fractional BEC, described by a spinor wave function with components (ϕ + , ϕ − ), which maintains the spin-orbit coupling (SOC) between the components [42].The respective system of coupled FGPEs in 2D is [37] with the kinetic-energy operator defined as per Eq. ( 14), cf.Eq. (19).Here, the attractive sign of the nonlinear interactions is assumed, to provide a possibility to create solitons.The coefficient of the self-attraction is scaled to be 1, while γ > 0 is the relative strength of the cross-attraction between the components, and λ is the strength of the linearly-mixing terms representing the SOC.By means of additional scaling, it is possible to set λ ≡ 1/2, while soliton states with real chemical potential µ are looked for as (cf.Eq. ( 15)), where complex stationary wave functions u ± satisfy equations Soliton solutions are characterized by the total norm, It is relevant to mention that the stationary system ( 22) can be derived from the respective Lagrangian: where * stand for the complex conjugate, and the integration with respect to p and q is reduced from the original domains (−∞, +∞) to (0, ∞), using the domains' symmetry.Although expression (24) seems cumbersome, it can be used to predict solitons solutions of Eq. ( 22) by means of the variational approximation (VA).To this end, the ansatz for the approximate solutions is adopted as with inverse width 1/ √ β and amplitudes A ± (cf. the known version of the VA for the usual 2D nonlinear Schrödinger equation [44,45]).The norm of ansatz (25) (25) represents the semi-vortex (SV) type, which implies that components u + and u − carry vorticities 0 and 1, respectively [43] (see an example shown by means of cross-sections of the two components in Fig. 1(b)).The SV structure is an exact feature of the 2D solitons shaped by SOC, which is not restricted by the applicability of VA.
The substitution of ansatz (25) in Lagrangian (24) produces the respective effective Lagrangian, Then, for given µ, SV's parameters are predicted by the Euler-Lagrange equations, ∂L eff /∂A ± = ∂L eff /∂β = 0. Figure 1(a) demonstrates that, as predicted by the VA and confirmed by numerical results, the fractional SOC system gives rise to stable SVs in the interval of LI values at norms below a critical value, N < N (α) the solitons are destroyed by the collapse, i.e., catastrophic self-compression of the field leading to the appearance of a singularity after a finite evolution time.The proximity of the VA-predicted and numerically found stability boundaries in Fig. 1(a) demonstrates that the VA, based on the simple ansatz (25), may produce reasonable results even for the complex system including the fractional diffraction (kinetic-energy operator), SOC, and the nonlinear interactions.
An exact property of the system is that N (SV) c (α = 1) = 0, i.e., no SVs exist at α ≤ 1.In the case of α − 1 → +0, the SV exists with an infinitesimal amplitude, hence it can be found as a numerical solution of the linearized version of Eq. ( 22) with α = 1.Cross-sections of components |u ± | of the latter solution are displayed in Fig. 1(b).Actually, these plots adequately represent the shape of the SV solitons in the general case, when the nonlinearity is present.
The SV solitons are stable for γ < 1 in Eqs.(20).For γ > 1, the SVs are unstable, while in this case there are stable solitons in the form of mixed modes, which mix zero-vorticity terms and ones with vorticities +1 and −1, see further details in Ref. [37].

III. THE EMULATION OF THE FRACTIONAL SCHR ÖDINGER EQUATIONS (FSES) IN OPTICS
A. Linear equations for the paraxial light propagation in optical systems emulating fractional diffraction Realizations of the fractional quantum systems were proposed in solid-state settings, such as Levy crystals [46] and exciton-polariton condensates in semiconductor microcavities [47].However, no experimental demonstration of the fractional quantum mechanics has been demonstrated thus far.
A more promising approach to the physical realization of FSEs in the form of classical equations is suggested by the commonly known similarity of the Schrödinger equation in quantum mechanics and the propagation equation for the amplitude of the classical optical field in the usual case of paraxial diffraction, with time t replaced, as the evolution variable, by the propagation distance, z [48].A scheme for the realization of this possibility was proposed by Longhi in 2015 [49], who considered the transverse light dynamics in optical cavities with the 4f (four-focal-lengths) structure.As shown in Fig. 2, the proposed setup incorporates two lenses and a phase mask, which is placed in the middle (Fourier) plane.The lenses perform the direct and inverse Fourier transforms of the light beam with respect to the single or two transverse coordinates, thus implementing Eq. ( 10) or (13) and the inverse form of these equations.The Fourier decomposition splits the beam into spectral components with transverse wavenumbers (p, q), the distance of each component from the optical axis being, roughly speaking, The action of the fractional diffraction in the Fourier space amounts, according to Eq. ( 9) or ( 14), to the local phase shift where Z is the propagation distance which accounts for the fractional diffraction.The corresponding differential phase shift of the Fourier components is imposed by the phase mask in Fig. 2, which is designed so that the local shift introduced my the mask at distance R from the axis emulates relation (29), with R and (p, q) related according to Eq. ( 28).In reality, the required mask, which is the central element of the setup, can be created as a computer-generated hologram [11].Thus, the light beam recombined by the right lens in Fig. 2 carries the phase structure corresponding to the action of the fractional diffraction.
In addition to that, the curved mirror at the left edge of the cavity in Fig. 2 introduces (prior to the action of the Fourier decomposition) a phase shift which represents the action of potential V (x, y) in Eq. (11).The layer of the gain medium placed next to the mirror in Fig. 2 may be used to amplify particular modes in the transverse structure of the light beam.
FIG. 2. The setup proposed in Ref. [49] for the emulation of the FSE in the optical cavity.The left and right lenses perform, respectively, the spectral (Fourier) decomposition and recombination of the light beam with respect to the transverse coordinates.The phase mask placed in the central Fourier plane introduces the differential phase shift which emulates the fractional diffraction,as per Eqs.( 28) and ( 29).The spherical mirror and gain medium at the left edge of the cavity introduce an effective potential.Reproduced with permission from S. Longhi, Fractional Schrödinger equation in optics, Opt.Lett.40, 1117-1120 (2015) [see Ref. [49]].Copyright 2015, Optica.
The setup which is outlined above provides a single-step transformation of the optical beam.The continuous FSE (11) is then introduced as an approximation for many cycles of circulation of light in the cavity, assuming that each cycle introduces a small phase shift (29).The 1D version of the FSE is obtained as the obvious reduction of the 2D scheme.
B. The emulation of the fractional diffraction in nonlinear optical systems: Fractional nonlinear Schrödinger equations (FNLSEs) in the spatial domain

The cubic nonlinearity in 1D
Once the FSE may be implemented as the propagation equation in optics under the action of the effective fractional diffraction, a natural possibility is to include the self-focusing nonlinearity of the optical material in the same system.The accordingly modified one-dimensional Eq. ( 8) is the fractional nonlinear Schrödinger equation (FNLSE), where g > 0 is the coefficient of the Kerr-nonlinearity coefficient, cf.FGPE (19).Then, steady-state solutions to Eq. ( 30), with real propagation constant k, are looked for as cf. Eq. ( 15), where real function U (x) satisfies the equation Localized solutions, i.e., fractional solitons, are characterized by their power (alias norm, cf.Eq. ( 23), The stability of these states against small perturbations is investigated by looking for solutions in the form of where a(x) and b * (x) represent the perturbation, and λ, which may be complex, is the instability growth rate.Substituting expression (34) in Eq. ( 30), one derives linearized equations for the small perturbations, The stationary solution (31) is stable if all eigenvalues produced by the numerical solution of Eq. ( 35) have Re(λ) = 0. Straightforward analysis of Eq. ( 30) reveals a scaling relation between the soliton's power and propagation constant: with some coefficient P 0 (α).Relation (36) satisfies the necessary stability condition, viz., the celebrated Vakhitov-Kolokolov (VK) criterion, dP/dk > 0 [50,51] at α > 1, suggesting that the corresponding soliton family may be stable.The case of α = 1, which corresponds to dP/dk = 0 in Eq. ( 36), implies the occurrence of the critical collapse, which makes all solitons unstable, similar to the family of Townes solitons produced by the 2D cubic nonlinear Schrödinger equation with the normal (non-fractional) diffraction (α = 2) [51][52][53].For α < 1, relation (36) contradicts the VK criterion, yielding dP/dk < 0, which implies strong instability of the solitons, driven by the supercritical collapse [51,52].
Similar to what is outlined above for the system of FGPEs (20), approximate fractal-soliton solutions of FNLSE (32) can be sought for by means of VA [21,54].To this end, the Lagrangian of Eq. ( 32) is used (cf.Eq. ( 24)), The variational ansatz for the solitons can be adopted as the usual Gaussian with width W and amplitude A (cf. Eq. ( 25)), whose power (33) is The substitution of ansatz (38) in Lagrangian (37) and integration yields the effective Lagrangian, in which the amplitude is eliminated in favor of the power by means of Eq. ( 39): cf. Eq. (26).Typical examples of the fractional solitons with the parameters predicted by the Euler-Lagrange equations, ∂L eff /∂P = ∂L eff /∂W = 0, and their comparison to numerical solutions obtained for the same values of α and k are displayed in Fig. 3. Further, the characteristic of the soliton family in the form of the P (k) dependence, as predicted by the VA and produced by the numerical solution, is plotted in Fig. 4.These results demonstrate that VA provides reasonable accuracy.
An important result of the VA is the prediction of the fixed power for the quasi-Townes solitons at α = 1, for which, as said above, the critical collapse takes place [54]:  The numerically found counterpart of the VA prediction ( 41) is Numerical investigation of the stability of the fractional-soliton family demonstrates that some of them may be weakly unstable at values of LI which are relatively far from α = 2.The instability does not destroy the solitons, leading to the onset of small-amplitude intrinsic oscillations in them [35].

Quadratic nonlinearity in 1D and 2D
In the case of quadratic (rather than cubic) self-focusing nonlinearity in the 1D space, the critical collapse takes place at α = 1/2, hence stable solitons may exist in the interval of 1/2 < α ≤ 2 [35], which is broader than its counterpart (27) in the case of the cubic self-focusing.In particular, the quadratic term may appear in Eq. ( 19) as −ε|ψ|ψ with ε > 0, representing a correction to the 1D FGPE induced by quantum fluctuations around the mean-field state [55].In optics, the well-known quadratic nonlinearity appears in the system of coupled FNLEs for the second-harmonicgeneration model: 2i where Ψ 1 and Ψ 2 are amplitudes of the fundamental and second harmonics, and real Q is the mismatch parameter, whose value may be fixed, by means of scaling, asQ = ±1 or 0. Soliton solutions to Eq. ( 43), with real propagation constants β and 2β of the fundamental-frequency and secondharmonic components, respectively, are looked for as where ψ 1,2 (x) are real localized profiles, see typical examples in Fig. 5.The existence and stability regions for solitons produced by system (43) in its parameter space were identified in Ref. [40].In particular, the stability condition completely coincides with the above-mentioned VK criterion, dP/dβ > 0, applied to the system's total power, P (β) = n=1,2 +∞ −∞ n 2 ψ 2 n (x)dx.In the 2D space, the cubic self-focusing nonlinearity gives rise to the critical collapse, as mentioned above, for α = 2, and to the supercritical collapse at all values of α < 2, therefore all the 2D fractional solitons are unstable.2D solitons, including ones with embedded vorticity, may be stabilized in the model including the cubic self-focusing and quintic defocusing.As shown in Ref. [26], families of stable solitons produced by the 2D FNLSE with the cubic-quintic nonlinearity and α < 2 are qualitatively similar to their counterparts which were studied in detail in the 2D equation with the same nonlinearity and normal (non-fractional) diffraction, corresponding to α = 2 [57].
On the other hand, the 2D version of the second-harmonic-generating system based on Eq. ( 43) initiates the collapse only at α ≤ 1, while 2D solitons may be stable in the above-mentioned intervals (27).However, such solutions have not been investigated, as yet.

Fast moving modes in the fractional medium: reduction to the non-fractional equation
The definition of the Riesz derivative, given by Eq. ( 9), breaks the commonly known Galilean invariance of Eq. ( 30) with the non-fractional diffraction (α = 2), in the absence of the external potential (V = 0).Nevertheless, the FNLSE (as well as the linear FSE) can be essentially simplified for fast moving modes (in the spatial domain the "moving" modes are actually ones tilted in the (x, z) plane) [12].To this end, one can consider solutions built as a product of a slowly varying amplitude ψ(x) and a rapidly oscillating continuous-waver (CW) carrier with large wavenumber P: The substitution of ansatz (45) in the nonlocally defined Riesz derivative in Eq. ( 9) leads to a quasi-local expression expanded in powers of small 1/P: Next, the substitution of expansion ( 46), truncated at a finite order, n = n max , in Eq. ( 30) produces the local nonlinear Schrödinger equation with higher-order diffraction terms corresponding to n ≥ 3 in Eq. (46).As concerns the convective term ∼ i∂ψ/∂x, which corresponds to n = 1 in Eq. ( 46), it may be eliminated by rewriting the resulting equation in the coordinate system moving (tilted) with the group velocity i.e., replacing coordinate x by x ≡ x − V gr z.Lastly, the term with n = 2 in expansion (46) gives rise to the usual second-order diffraction term, with coefficient . Thus, the resulting equation, truncated at n max = 2, takes the form of the usual (non-fractional) nonlinear Schrödinger equation, which gives rise to the commonly known solutions, such as solitons.
C. Systems of fractional nonlinear Schrödinger equations with the cubic nonlinearity

Fractional domain walls
Extension of the fractional model for copropagating optical waves with orthogonal polarizations or different carrier wavelength is a natural possibility, which was addressed in Ref. [36] for the case of the self-defocusing cubic nonlinearity.The corresponding 1D system of coupled FNLSEs for wave amplitudes Ψ and Φ is where β > 0 is the relative strength of the cross-phase-modulation (XPM) interaction, while the coefficient of the self-phase modulation is scaled to be 1.In the case of orthogonal polarizations, natural XPM values are β = 2/3 or 2, for the linear and circular polarizations, respectively.The copropagation of optical waves carried by different wavelengths is governed by the same system (49) with β = 2 [48].Other positive values of β are possible in photonic crystals [48].
In addition to XPM, the two components may be coupled in Eq. ( 49) by linear terms with real coefficient λ.The linear coupling accounts for the twist of the waveguide in the case of the copropagation of linear polarizations, or elliptic deformation of the waveguide in the case of circular polarizations [48].
As usual, stationary states with a real propagation constant k are looked for as {Ψ, Φ} = exp (ikz) {U (x), V (x)}, where real real functions U (x) and V (x) are solutions of the system Natural patterns supported by Eq. ( 50) with k < 0 are optical domain walls (DWs) [58][59][60], which link different CW states at x → ±∞, that are mirror images to each other:  [36]] with the permission of AIP Publishing.
with equal total power densities: In particular, in the absence of the linear coupling, λ = 0, U (x) and V (x) vanish at x → +∞ and −∞, respectively.The asymmetric CW states (51) exist under the condition of the immiscibility [61] of the two components: In the case of β = 3 (which plays a special role in various systems, making it possible to find exact DW solutions [60,62]), substitution {U (x), V (x)} = (1/2) √ λ − k ∓ W (x) reduces two equations ( 50) to a single one, In this case, the boundary conditions (51) which determined the DW solutions reduce to Note that it follows from Eq. ( 52) with β = 3 that √ −k − λ is real, hence the boundary conditions (54) make sense.The DW states produced by Eq. ( 49) are stable (in particular, the CW background states at x → ±∞ are not subject to modulational instability).A set of numerically obtained DW solutions are displayed in Fig. 6.
A more general bimodal system, with different values of LI in its coupled components, α 1 and α 2 , was considered too [36].It also gives rise to families of stable DW patterns, which, unlike those displayed in Fig. 6, are spatially asymmetric with respect to the midpoint.

Couplers and spontaneous symmetry breaking
Another model which naturally gives rise to a system of coupled FNLSEs with the fractional diffraction and cubic self-focusing acting in each equation, and linear coupling between the equations, introduces a dual-core optical waveguide [34,38]: where the coefficient of the linear coupling is scaled to be 1 (cf.coefficient λ in Eq. ( 49)).As well as the single FNLSE (30), system ( 55) is meaningful in the interval of LI index given by Eq. ( 27), when the collapse does not take place.Steady-state solutions to Eq. ( 55) with real propagation constant k > 0 are looked for Ψ 1,2 = exp (ikz) U 1,2 (x), where real functions U 1,2 (x) satisfy equations Two-component solitons produced by Eq. ( 55) are characterized by the total power, In the case of the normal (non-fractional) diffraction, α = 2, a fundamental property of the two-component solitons is spontaneous symmetry breaking (SSB), which destabilizes the obvious symmetric solitons, with U 1 (x) = U 2 (x), at a critical value of the soliton's power, P = P cr , and replaces the unstable symmetric modes by asymmetric ones, with U 1 (x) ̸ = U 2 (x).The asymmetry of the solitons is naturally quantified by the normalized difference in powers of their components, i.e., The SSB phenomenology for solitons was studied in detail theoretically [63]- [68], and was recently demonstrated experimentally in dual-core optical fibers [69] with the second-order (non-fractional) GVD.A characteristic peculiarity of the SSB in this system is that the transition from the symmetric solitons to asymmetric ones is performed through a subcritical bifurcation [70] (similar to a symmetry-breaking phase transition of the first kind), which means that a stable asymmetric soliton solution emerges at a value of P which is slightly smaller than P cr , i.e., in the subcritical range.In this case, the branch of the asymmetric soliton states, appearing at P = P cr , originally moves in the backward direction, being unstable (therefore, another name of this bifurcation is the backward one).Then, this branch turns forward, getting stable, at the turning point (see, e.g., the left panel of Fig. 7).
The system of equations (56), with the fractional Riesz derivative, defined as per Eq. ( 9), can be derived from the Lagrangian, where the Hamiltonian is cf. Eq. ( 37).Using the Lagrangian, VA can be applied to the system.An appropriate tractable ansatz, with width W , power P , and angle χ accounting for the norm distribution between the components, is In terms of this ansatz, the asymmetry parameter ( 58) is The calculation of Lagrangian (59) with ansatz (61) yields (cf.Eq. ( 40)), where Γ and ζ are the Gamma-and zeta-functions.The effective Lagrangian (63) gives rise to the Euler-Lagrange equations, ∂L eff /∂W = ∂L eff /∂ (sin(2χ)) = 0, which amount to relation W = P S/6 and an equation for parameter S (which is defined in Eq. ( 62)), The SSB point, as which the family of the asymmetric solitons branches off from the symmetric one, corresponds to setting S = 1 (which represents the symmetric soliton) in Eq. ( 64), the result being In the case of α = 2, Eq. ( 65) reproduces the known result produced by the VA for the non-fractional coupler, viz., [67,68].Note that the SSB point for the usual coupler is known in an exact form, which can be found without the use of VA [63], i.e., the relative error of the VA prediction is ≈ 5.7%.In the opposite limit of α − 1 → +0, Eq. ( 65) gives the respective VA result as (P SSB ) (α → 1) = (12/π) ln 2 ≈ 2.648.Eventually, the VA predicts the basic characteristic of the system, viz., the dependence of the asymmetry parameter (58) on the power, Θ(P ), in an implicit form [38]: In the limit of α → 1, this relation yields The VA predictions following from Eqs. ( 67) and ( 65) for different values of LI α are plotted in Fig. 7, along with results produced by the numerical solution of Eq. (56).It is observed that the subcritical character of the SSB bifurcation becomes more pronounced with the decrease of α.In the limit of α = 1, the dependence is available solely in the VA form, given by Eq. (68), while the numerical solutions were obtained only for α ≥ 1.2 [38].The VA-predicted Θ(P ) dependence for α = 1, given by Eq. ( 68), takes the extreme subcritical form (defined as per Ref. ( [39]), in which the branch of the asymmetric solitons, going backward from the SSB point, never turns forward, always remaining unstable).The dependence (68) terminates at point Θ = 1, i.e., at P = (6/π) ln 2 ≈ 1.324 (see the left panel in Fig. 7), as definition (58) does not admit values |Θ| > 1.
In addition to that, Fig. 8 displays a set of numerically generated dependences P (k) for different fixed values of LI.The SSB points are those separating stable and unstable segments of the symmetric-soliton branches, from which ones representing unstable asymmetric solitons stem.The stability boundaries, observed on the asymmetric branches for α = 2.0, 1.8, 1.6, and 1.4, correspond to the above-mentioned turning points, while the branches for α = 1.2 and 1.1 do not reach those points, therefore they seem completely unstable in Fig. 8.
The analysis of the SSB phenomenology was also developed for moving (tilted) two-component solitons (recall it is a nontrivial extension of the analysis because the fractional diffraction breaks the system's Galilean invariance), demonstrating that the symmetry-breaking bifurcation keeps its subcritical character, and values of the power at the bifurcation point become lower than for the quiescent solitons [38].Collisions between moving solitons were studied too, with a conclusion that the collisions with small relative velocities lead to elastic rebound, which is followed by strongly inelastic symmetry-breaking interactions at intermediate values of the velocities, and, eventually, by restoration of the elasticity in collisions between fast solitons [38].
The analysis of SSB in solitons was also developed for the single-component FNLSE (30) with a symmetric doublewell potential V (x) [33].On the contrary to what is shown here, in that case, the SSB bifurcation is of the supercritical (forward ) type (i.e., it represents a symmetry-breaking phase transition of the second type).Furthermore, the same FNLSE with the self-defocusing nonlinearity (g < 0 in Eq. ( 30)) does not break the symmetry of the ground state trapped by the double-well potential, but the increase of |g| leads to spontaneous breaking of the antisymmetry of the first excited state trapped by the same potential.57), as predicted by the VA through Eq. ( 67), and as produced by the numerical solution of Eq. ( 56), for different values of LI α.Right: the value of the soliton's power at the SSB bifurcation point vs. the LI, as predicted by the VA, and as found from the numerical solution.
The red circle at the end of the numerical curve represents the known exact value (66) for α = 2 (the usual system with the non-fractional diffraction).Reproduced from D. V. Strunin and B. A. Malomed, Symmetry-breaking transitions in quiescent and moving solitons in fractional couplers, Phys.Rev. E 107, 064203 (2023) [see Ref. [38]] with the permission of AIP Publishing.Thus far, no experimental realization of the effective fractional diffraction for light beams in the spatial domain has been reported.Another option for the implementation of the concept of the fractional propagation in optics is to resort to the transmission of light pulses in the temporal domain, i.e., in fiber-based setups, with an effective fractional GVD.This option has been realized recently [11], thus providing the first experimental implementation of a fractional FIG. 9. (a) The setup elaborated in Ref. [11] for the realization of the fractional GVD in the fiber-laser cavity.The hologram installed in the initial section is employed to shape the input pulse.In the propagation section, the second hologram plays the role of the phase mask (cf.medium in any physical setting. Confining the consideration to the linear propagation regime, the effective FSE for the optical amplitude Ψ (τ, z) in the temporal domain is written as where the evolutional variable is the propagation distance z, and the effective coordinate is the reduced time [48], τ = t−z/V gr , where V gr is the group velocity of the carrier wave.In this equation, the fractional GVD is determined by LI α (cf.Eq. ( 8)), D is the corresponding coefficient, and the usual derivatives represent the second-and higher-order GVD, for k = 2 and k ≥ 3, respectively, which act in the setup along with the fractional GVD.The solution of Eq. ( 69) can be easily written for the Fourier transform of the optical field, (cf.Eq. ( 10)), as where Ψinput (ω) is the Fourier transform of initial condition, i.e., the optical signal coupled into the fiber at z = 0.The fiber-cavity setup which was used for the experimental emulation of the fractional GVD is displayed in Fig. 9.It incorporates two holograms, with the left one shaping the input pulse.The second hologram, installed at the central position, was designed as the element similar to the phase mask in the setup presented above in Fig. 2 (it is displayed in the rotated form, to make its internal structure visible).It imposes the differential phase shift onto spectral components of the light signal, simulating the action of the fractional GVD as per Eq.(71) (cf.Eq. ( 29)).
Experimental findings were compared with the theoretical result provided by Eq. ( 71), keeping only the secondorder GVD ∼ β, in addition to the fractional GVD [11].Typical values of the dispersion coefficients in Eq. ( 71) where β 2 < 0 implies that it represents the anomalous second-order GVD [48], although the case of the normal GVD, corresponding to β 2 > 0, was considered too.The results were summarized as functions of two control parameters, namely, LI α and effective GVD length L GVD (which corresponds to z in Eq. ( 71)).The plane of these parameters is schematically defined in the inset to Fig. 9, where α varies from 0 to 2, while positive and negative values of L GVD correspond to the anomalous and normal second-order GVD.Basic experimental and theoretical results are reported in Fig. 10 for four characteristic cases, with parameters (L GVD , α) = (5, 1.25), (5, 0.25), (−5, 0.25), (−5, 1.25), that pertain, respectively, to quadrants Q1 -Q4 in Fig. 9(b).In Q1, the setup demonstrates competition between the splitting effect of the fractional GVD and the confining (anti-splitting) action of the anomalous second-order GVD.The competition maintains balance over a short propagation distance, ≃ 40 m.Then, the splitting effect becomes the dominant one, leading to fission of the pulse in two secondary ones.In Q2, small α = 0.25 implies that the fractional GVD is weak.Then, the dominant anomalous second-order GVD drives fragmentation of the input into a multi-jet pattern, which, however, stays confined, following the confining effect of the same GVD term observed at the initial stage of the evolution in Q1.In panel Q3, the fractional GVD remains weak, with the same value of LI, α = 0.25, as in Q2, while the second-order GVD flips its sign from anomalous to normal, leading to the fast split of the pulse in two.Finally, in quadrant Q4 the same normal second-order GVD as in Q3 cooperates with the strong fractional-GVD term (α = 1.25), producing extremely fast fission of the pulse in two fragments (note the difference by an order of magnitude between scales on the vertical axes in panels Q3 and Q4).Another difference between the relatively gentle splitting in Q3, driven by the normal second-order GVD, and more violent splitting in Q4, driven jointly by the second-order and fractional GVD terms, is that, in the latter case, the emerging fragments are loosely bound pulses, unlike tightly bound ones in Q3.The dynamics displayed in Fig. 10 is reversible, therefore examples of the splitting, such as ones exhibited in panels Q1 and Q3, imply the possibility of fusion of colliding pulses into a single one, if the propagation distance z varies in the opposite direction.
For the comparison's sake, panels B2 in Fig. 10 demonstrate that, unlike the setup including the fractional GVD, the usual one, with L GVD = 0 and α = 2 (see area B2 in the inset in Fig. 9(b)), very quickly leads to full destruction of the input pulse (note that largest propagation distance in panels B2 is 0.02 km).Thus, the fractional GVD makes the outcome of the light propagation in the fiber cavity much more nontrivial.
A remarkable fact is close proximity between the experimental findings and the corresponding theoretical results, which are presented, respectively, in rows (b) and (a) of Fig. 10.In addition to that, Fig. 11 reports a set of theoretical results for the same dispersion parameters as in Eq. ( 72), but with another value of the LI, α = 1.In particular, Fig. 11(a) demonstrates that, under the action of fractional-only GVD (L GVD = 0), the evolution of the input is generally similar to that in the case when the usual second-order GVD is present too, but the fractional term is the dominant one (case Q1 in Fig. 10): the input pulse quickly splits in a pair of secondary ones.Such a setting with the purely fractional GVD was not realized experimentally, as the usual GVD cannot be eliminated in the real setup.Further, if strong anomalous or normal second-order GVD is included, corresponding to L GVD = +10 and −10 in Figs.11(b) and (c), respectively, the outcomes are also akin to what is shown for similar cases in panels Q2 and Q4 of Fig. 10: the fragmentation of the input into a non-expanding multi-jet pattern in the former case, or violent fission of the input in two loosely bound pulses in the latter case.The experiments reported in Ref. [11] were also performed for the input pulse with the phase including a cubic term, ∼ τ 3 .It is well known that, in the framework of the usual linear Schrödinger equation with the second-order GVD, this term initiated the propagation of a self-accelerating self-bending Airy wave [71,72].In the interval of LI values 1 ≤ α ≤ 1.80, the results reported in Ref. [11] (not shown here in detail) demonstrate splitting of the input into a self-bending quasi-Airy wave and an additional one propagating along a straight trajectory.
All the experimental results demonstrating the implementation of the effective fractional GVD, summarized above and reported in detail in Ref. [11], have been obtained in the linear regime.The interplay of the strongly fractional GVD, with LI values α = 1 and α = 0.2, and the fiber's self-focusing nonlinearity is a subject of a very recent experimental work [82].In particular, it demonstrates spectral bifurcations of ultrashort optical pulses in the fiberlaser cavity, in the form of spontaneous transitions between the pulses with single-and multi-lobe structures.

V. CONCLUSION
This aim of this paper is to present a concise summary of dynamical models of 1D and 2D media with fractional diffraction or dispersion, both linear and nonlinear ones.Two well-substantiated physical models of this type are fractional quantum mechanics for particles which, at the classical level, move by random Lévy flights [13][14][15], and the emulation of the effective fractional diffraction in optical cavities [49].These models include the Riesz fractional derivatives [10] or the respective 2D fractional Laplacian, which are defined as per Eqs.( 9) and ( 14).The Riesz derivative with LI (Lévy index) α amounts to the multiplication of the Fourier transform of the underlying wave field, Ψ(p), by |p| α .Equations ( 9) and ( 14) demonstrate that the Riesz derivatives and their 2D counterparts are represented, in the explicit form, not by differential operators, but rather by nonlocal integral ones.
The paper presents basic types of solitons produced by the nonlinear fractional models in a brief form, without an objective to produce a systematic review of the solitons in fractional systems.The absolute majority of theoretical results for the fractional solitons have been obtained by means of numerical methods.Nevertheless, the present review includes some quasi-analytical results based on the VA (variational approximation).Actually, the comparison with numerical findings demonstrates that the VA works surprisingly well too in these relatively complex nonlinear nonlocal models.
In addition to the theoretical results, the article includes a summary of the recently reported first experimental realization of the fractional optics.It is based on the fiber-cavity setup, in which the effective fractional GVD (groupvelocity dispersion) is implemented by means of a specially designed phase plate (computer-generated hologram) [11].
There are essential topics dealing with the dynamics of fractional media which are not considered in this paper, which is not designed as a comprehensive review.One of such topics is the dynamics of fractional discrete systems.[73][74][75].Others, which have attracted much interest recently (and may be relevant subjects for a more extensive review), are dissipative solitons and self-trapped vortices in models based on fractional complex Ginzburg-Landau equations [54,[76][77][78], as well as solitons in fractional parity-time-symmetric systems [79][80][81].
There are many possibilities for further theoretical and, especially, experimental developments in this area.In particular, on the theoretical side, a challenging objective (which is mentioned in this paper) is to derive equations of the FGPE (fractional Gross-Pitaevskii equation) type for a BEC (Bose-Einstein condensate) of quantum particles which are individually governed by the fractional Schrödinger equation.On the experimental side, obvious objectives are to realize the fractional diffraction in spatial-domain optics, and to create the predicted fractional solitons in such settings.

FIG. 1 .
FIG. 1.(a) The 2D SV solitons exist and are stable at norms below the critical value, N (SV) c , which is plotted vs. the half-LI, α/2, for λ = 0.5 and γ = 0 in Eq. (4).The collapse takes place at N > N (SV) c .The numerical result and its VA-predicted counterpart are shown by the dashed line and chain of rhombuses, respectively.(b) Cross sections of profiles of |ϕ+| and |ϕ−| in the SV state (continuous and dashed lines, respectively) produced by the numerical solution of the linearized version of Eq. (22) with α = 1 (see the text).Reproduced from H. Sakaguchi and B. A. Malomed, One-and two-dimensional solitons in spin-orbit-coupled Bose-Einstein condensates with fractional kinetic energy, J. Phys.B: At.Mol.Opt.Phys.55, 155301 (2022) [see Ref. [37]], with the permission of IOP Publishing.

FIG. 6 .
FIG. 6.A set of stable DW patterns with k = −1, produced by the numerical solution of Eq. (50) in the absence of the linear coupling (λ = 0), with LI α = 1, and different values of the XPM coefficient β indicated in the figure.Reproduced from S. Kumar, P. Li, and B. A. Malomed, Domain walls in fractional media, Phys.Rev. E 106, 054207 (2022) [see Ref. [36]] with the permission of AIP Publishing.

FIG. 7 .
FIG. 7.Left: A set of dependences of the soliton's asymmetry parameter (58) on total power(57), as predicted by the VA through Eq. (67), and as produced by the numerical solution of Eq. (56), for different values of LI α.Right: the value of the soliton's power at the SSB bifurcation point vs. the LI, as predicted by the VA, and as found from the numerical solution.The red circle at the end of the numerical curve represents the known exact value (66) for α = 2 (the usual system with the non-fractional diffraction).Reproduced from D. V. Strunin and B. A. Malomed, Symmetry-breaking transitions in quiescent and moving solitons in fractional couplers, Phys.Rev. E 107, 064203 (2023) [see Ref.[38]] with the permission of AIP Publishing.

3 FIG. 8 .
FIG.8.Dependences between the total power, P , and the propagation constant, k, for families of symmetric and asymmetric solitons produced by the numerical solution of Eqs.(55) and (56), for different fixed values of LI α.Solid and dashed segments of the curves represent, respectively, stable and unstable parts of the soliton families.Reproduced from D. V. Strunin and B. A. Malomed, Symmetry-breaking transitions in quiescent and moving solitons in fractional couplers, Phys.Rev. E 107, 064203 (2023) [see Ref.[38]] with the permission of AIP Publishing.

Fig. 2 )
FIG. 9. (a)The setup elaborated in Ref.[11] for the realization of the fractional GVD in the fiber-laser cavity.The hologram installed in the initial section is employed to shape the input pulse.In the propagation section, the second hologram plays the role of the phase mask (cf.Fig. 2), which imposes the differential phase shift onto spectral components of the light signal, to emulate the effect of the fractional GVD.This hologram is shown in the rotated position, to display its structure.The initial temporal profile of the pulse belonging to quadrant Q3 in (b), and its evolution in the propagation section are shown in the top of the panel.(b) Four different quadrants in the parameter plane (LGVD, α).Q1 and Q2 correspond to the cases of LGVD > 0 with 1 ≤ α ≤ 2 and 0 ≤ α ≤ 1, respectively.Q3 and Q4 are defined similarly for LGVD < 0. Areas B1 and B2 designate narrow strips with α close to 0 or 2, respectively.Reproduced with permission from S. Liu, Y. Zhang, B. A. Malomed, and E. Karimi, Experimental realisations of the fractional Schrödinger equation in the temporal domain, Nature Comm.14, 222 (2023) [see Ref. [11]].Copyright 2023, Springer-Nature.

FIG. 10 .
FIG. 10.Different outcomes of the propagation of input light pulses in parameter regions Q1 -Q4 and B2, defined as per Fig. 9(b), are displayed in the plane of the temporal coordinate τ and propagation distance z, see Eq. (69).The color code represents local intensities of the pulses.Exact values of the parameters corresponding to Q1 -Q4 are fixed in Eq. (73), while B2 corresponds to LGVD = 0 and α = 2.In all the cases, the dispersion parameters are fixed as per Eq.(72).Rows of panels (a) and (b) represent, severally, theoretical and experimental results.Note the great difference in the scales of the vertical axes (z) in panels Q2, Q3 (zmax = 1 km), Q1, Q4 (zmax = 0.1 km) and B2 (zmax = 0.02 km).Reproduced with permission from S. Liu, Y. Zhang, B. A. Malomed, and E. Karimi, Experimental realisations of the fractional Schrödinger equation in the temporal domain, Nature Comm.14, 222 (2023) [see Ref. [11]].Copyright 2023, Springer-Nature.

FIG. 11 .
FIG. 11.Different outcomes of the propagation of the input pulses, as predicted by simulations of the model based on Eq. (71) with the GVD parameters fixed as per Eq.(72), LI α = 1, and values of LGVD indicated in panels (a,b,c), are displayed in the plane of (τ, z).The color code represents local intensities.Reproduced with permission from S. Liu, Y. Zhang, B. A. Malomed, and E. Karimi, Experimental realisations of the fractional Schrödinger equation in the temporal domain, Nature Comm.14, 222 (2023) [see Ref. [11]].Copyright 2023, Springer-Nature.