Physica D 462 (2024) 134157

Contents lists available at ScienceDirect

Physica D
journal homepage: www.elsevier.com/locate/physd

Formation and propagation of fundamental and vortex soliton families in 1D

and 2D dual-Lévy-index fractional nonlinear Schrödinger equations with
cubic–quintic nonlinearity
Ming Zhong a,b , Yong Chen c , Zhenya Yan a,b ,∗, Boris A. Malomed d,e
a
KLMM, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China
b School of Mathematical Sciences, University of Chinese Academy of sciences, Beijing 100049, China
c School of Mathematics and Statistics, Jiangsu Normal University, Xuzhou 221116, China
d Department of Physical Electronics, School of Electrical Engineering, Faculty of Engineering, Tel Aviv University, Tel Aviv 69978, Israel
e
Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica, Chile

ARTICLE INFO ABSTRACT

Keywords: We address effects of interplay of two different fractional-diffraction terms, corresponding to different Lé vy
1D and 2D dual-Lévy-index fractional indices (LIs), in the framework of nonlinear Schrödinger equation (NLSE) with cubic–quintic nonlinearity (dual-
nonlinear Schrödinger equations LI fractional NLSE ), in one- and two-dimensional (1D and 2D) settings. The critical (in 1D) and supercritical (in
Cubic–quintic nonlinearity
2D) wave collapses are suppressed in the presence of a defocusing quintic term, making it possible to produce
Variational approximation
stable localized modes, including fundamental and vortical ones. In 2D, families of fundamental and vortex
Splitting instability
Anisotropic fractionality
solitons, with topological charges 𝑆 = 0, 1, 2, 3, are produced in a numerical form and by dint of the variational
approximation (VA). In particular, the threshold power necessary for the existence of 2D fundamental solitons
is predicted by the VA, being close to the numerical results. In the case of fractional-diffraction terms with
unequal LIs acting in two transverse directions, anisotropic 2D fundamental and vortex solitons are constructed.
Stability of the solitons is investigated for small perturbations governed by the linearized equations, and results
are corroborated by direct simulations of the perturbed evolution. Those vortex solitons which are unstable are
split by azimuthal perturbations into fragments, whose number is determined by the shape of the perturbation
eigenmode with the largest growth rate. These results extend the concept of the 2D fractional diffraction and
solitons to media with the anisotropic fractionality.

1. Introduction Eqs. (9) and (13)], and 𝑈 (𝐫, 𝑡) with 𝐫 = (𝑥1 , 𝑥2 , … , 𝑥𝑑 ) is an external
potential. Relying upon the commonly known similarity between the
Fractional quantum mechanics was first proposed by Laskin [1] in quantum-mechanical Schrödinger equation and the classical equation
2000 as an extension of the canonical quantum mechanics, in which for the paraxial propagation of light, with time 𝑡 replaced by propaga-
the Feynman-path integration is performed over random Lévy flights of tion distance 𝑧, Longhi [8] had proposed to implement the fractional
classical particles [2,3], instead of the usual Brownian trajectories. The
Schrödinger equation in optics, using the direct and inverse Fourier
respective fractional Schrödinger equation for wave function 𝜓 [4,5]
transforms of the laser beam in the transverse direction, with a specially
includes the fractional kinetic-energy (diffraction) term characterized
designed phase mask inserted in between, which imparts phase shifts to
by the Lévy index (LI) 𝛼 ∈ (1, 2] [6]
the Fourier components, emulating the action of the fractional diffrac-
𝜕𝜓 [ ( 2 2 )𝛼∕2 ]
𝑖ℏ = 𝐷𝛼 −ℏ ∇ + 𝑈 (𝐫, 𝑡) 𝜓. (1) tion [8–10]. Schemes for the implementation of the fractional quantum
𝜕𝑡
mechanics in condensed-matter settings were proposed too [11,12].
Here 𝐷𝛼 > 0 is a diffraction coefficient, and operator (−ℏ2 ∇2 )𝛼∕2 ,
Recently, the first experimental realization of the fractional group-
with the reduced Planck’s constant ℏ and Laplacian ∇2 = 𝜕𝑥2 + 𝜕𝑥2 +
1 2 velocity dispersion (GVD) in a fiber-laser cavity was reported [10].
⋯ + 𝜕𝑥2 , is defined as the Riesz fractional derivatives [7], carrying
𝑑 The scheme makes use of two pulse-shaping holograms, one of which
over into the usual second derivative in the case of 𝛼 = 2 [the
explicit definition of the fractional derivatives is presented below in produces an appropriate input, and the other one emulates the phase

∗ Corresponding author at: KLMM, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China.
E-mail address: zyyan@mmrc.iss.ac.cn (Z. Yan).

https://doi.org/10.1016/j.physd.2024.134157
Received 6 December 2023; Received in revised form 28 February 2024; Accepted 2 April 2024
Available online 6 April 2024
0167-2789/© 2024 Elsevier B.V. All rights reserved.
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

shifts corresponding to the action of the fractional dispersion in the and 𝑦 (actually, both structures may be integrated into a single phase
temporal domain. mask fabricated as an appropriate hologram). As for the 1D setting
The implementation of fractional diffraction in a nonlinear optical modeled by Eq. (4), a possibility to experimentally implement the dual-
medium gives rise to the fractional nonlinear Schrödinger equation LI system is suggested by the above-mentioned experimental work [10].
(FNLSE) for the slowly varying envelope 𝜓 of the optical field, see Indeed, aiming to implement the action of different fractional-GVD
reviews in Refs. [13,14], and references therein: terms, one may design a fiber cavity with two inserted holograms,
𝜕𝜓 ( )𝛼∕2 ( ) each representing its LI value. Furthermore, structures emulating both
𝑖 = 𝐷𝛼 −∇2 𝜓 + 𝑈 (𝐫, 𝑧)𝜓 + 𝑓 𝐫, |𝜓|2 𝜓, 𝛼 ∈ (1, 2], (2)
𝜕𝑧 fractional terms can be inscribed on the same hologram, which was
( )
where 𝑓 𝐫, |𝜓|2 𝜓 represents the nonlinearity, which may include actually implemented in Ref. [10], where the single hologram was
cubic and quintic terms, −𝑔|𝜓|2 𝜓 and 𝛾|𝜓|4 𝜓, with 𝑔 > 0 (< 0) used to control both the fractional GVD, with the set of LI values
and 𝛾 < 0 (> 0) corresponding to the self-focusing and defocusing, 𝛼 = 0.05, 0.25, 0.4, 0.6, 0.8, 1.0, 1.25, 1.5, and 1.8, and regular GVD, cor-
respectively. Many varieties of solitons have been predicted in such a responding to 𝛼 = 2. These options include the pair of the LI values
fractional model, including spatiotemporal “accessible solitons" [15], 𝛼1 = 2, 𝛼2 = 1.5, which plays the main role in the present work.
gap [16–18] and multi-pole modes [19,20], and solitary vortices in The objective of this work is to study families of fundamental and
multi-dimensional settings [21–23]. In particular, a noteworthy finding vortex solitons as solutions of 1D and 2D FNLSEs [in particular, Eq. (4)
is that, when considering the symmetry breaking phenomenon in the in 2D] with the dual-LI diffraction terms. These solutions are produced
FNLSE, it may give rise to “ghost" modes [24,25], in the form of in the numerical form and by means of a variational approximation
asymmetric solitons with a complex propagation constant [26–29]. (VA), demonstrating agreements between them, except for the case of
Solitons in the two-component system with the quadratic nonlinear- the near-threshold propagation. Further, we demonstrate that quasi-
ity, which models the second-harmonic generation in the presence of Townes solitons, which are generated by the 1D NLSE with the normal
the fractional diffraction, were predicted too [30]. Domain walls in diffraction (LI = 2) and quintic-only self-focusing nonlinearity [49–
fractional media were also considered recently [31].
51], are stabilized under the action of the second diffraction term.
It is well known that two-dimensional (2D) and three-dimensional
The 2D model generates isotropic and anisotropic fundamental solitons,
(3D) NLSEs with the normal (non-fractional) diffraction and cubic self-
as well as solitary vortices with integer topological charges (winding
focusing give rise to the critical and supercritical wave collapse, which
numbers) 1 ≤ 𝑆 ≤ 3. Stability areas for the fundamental and vortex
occurs, respectively, above a finite or zero threshold value of the input’s
solitons are established by dint of the linear-stability analysis and direct
norm, i.e., 2D power or 3D energy, in terms of optics [32,33]. The
collapse destabilizes solitons produced by these equations, including simulations.
the degenerate family of 2D Townes solitons (TSs), which exist at the The main findings reported in this work are summarized as follows:
single value of the power [34,35]. A common way to arrest the collapse
and stabilize the solitons is to add the defocusing quintic term, i.e., to • We present numerically found and VA-predicted families of soli-
consider the model based in the cubic–quintic (CQ) NLSE [36–42], ton solutions of the 1D dual-LI cubic–quintic FNLSE (6).
• The effect of quintic coefficient 𝛾 is analyzed for the solitons of
𝜕𝜓
𝑖 = −∇2 𝜓 − 𝑔|𝜓|2 𝜓 + 𝛾|𝜓|4 𝜓, (3) the 1D dual-LI cubic–quintic FNLSE (6).
𝜕𝑧
• We present numerically found and VA-predicted families of fun-
with 𝑔 > 0 and 𝛾 > 0. Eq. (3) produces families of fundamental and
damental and vortex solutions, with topological charges 𝑆 =
vortex solitons, which are partly stable. The 2D version of Eq. (3)
0, 1, 2, 3, of two types of 2D dual-LI cubic–quintic FNLSEs, viz., the
applies to optics, considering the light propagation in a bulk waveguide
ones given by Eqs. (11) and (14).
filled by a dielectric material whose intrinsic nonlinearity is modeled by
• We find the stable anisotropic 2D fundamental and vortex solitons
the CQ nonlinearity, such as liquid carbon disulfide [43] or a colloidal
suspension of nanoparticles [44]. in the case of fractional-diffraction terms with unequal LIs acting
It is also relevant to mention that combinations of fourth-order in the transverse directions. To the best our knowledge, these
or higher-order dispersion terms with the normal second-order one, found 2D soliton phenomena did not appear in the standard 2D
which occur in various optical systems [45,46], give rise to novel NLS equation with cubic and/or quintic nonlinearity.
species of temporal solitons. In particular, bright solitons arising from • We demonstrate that the critical and supercritical wave collapses
the interplay of the negative fourth-order dispersion and self-focusing are suppressed in the presence of the self-defocusing quintic term.
were observed experimentally [47]. Following this line of studies, it • The stability of the fundamental and vortex solutions is explored
is relevant to consider the formation and dynamics of solitons in the in the framework of the 1D and 2D dual-LI cubic–quintic FNLSEs.
framework of a dual-LI FNLSE, featuring a combination of fractional-
diffraction terms with two different values of LI, 𝛼1 and 𝛼2 (𝛼1,2 ∈ The rest of this paper is arranged as follows. We introduce the dual-
(1, 2]), i.e., LI CQ-FNLSE models in the detailed form, along with some methods
𝜕𝜓 [ ( )𝛼 ∕2 ( )𝛼 ∕2 ] ( ) necessary to work with them, in Section 2. 1D fundamental solitons
𝑖 = 𝐷𝛼1 −∇2 1 + 𝐷𝛼2 −∇2 2 𝜓 + 𝑓 𝐫, |𝜓|2 𝜓 (4) are produced by means of the VA and numerical methods in Section 3.
𝜕𝑧
[cf. Eq. (2)], where 𝐷𝛼𝑗 (𝑗 = 1, 2) are the respective real-valued We specifically explore effects of variation of the system’s parameters
diffraction coefficients. In this connection, it is relevant to mention that on properties of the solitons. In Section 4, we address the fundamental
spatiotemporal FNLSE of a mixed type, with the fractional diffraction and vortical solitons in the 2D settings, both isotropic and anisotropic
and ordinary (non-fractional) GVD acting in the spatial and temporal ones. The paper is concluded by Section 5.
directions, respectively, was introduced recently, and its modulational
instability and collapse dynamics were investigated analytically and 2. 1D and 2D dual-Lévy-index nonlinear models
numerically [48].
A natural possibility to introduce fractional-diffraction terms with
2.1. The basic equations
different LIs is suggested by the 2D FNLSE with unequal LIs acting along
two spatial directions, 𝑥 and 𝑦. In terms of the possible experimental
realization [8,10], this scheme implies the use of two phase masks We start the analysis by considering the propagation of laser beams
with mutually perpendicular quasi-1D structures, which impart the along the 𝑧 axis in a waveguide with transverse coordinate 𝑥 under the
differential phase shifts emulating the fractional diffraction along 𝑥 action of two fractional-diffraction terms and CQ nonlinearity, modeled

2
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

( )𝛼 ∕2
by the following CQ-FNLSE for complex amplitude 𝛹 (𝜉, 𝜁) of the light where operator −∇2 𝑗 is defined with the help of the 2D Fourier
field, written in the scaled form [22,36,37,52]: transform:
( 2 )𝛼∕2 [ ]
⎧ ( )𝛼1 ∕2 ⎫ 𝜓(𝑥, 𝑦) = 2−1 (𝑝2 + 𝑞 2 )𝛼∕2 2 (𝜓(𝑥, 𝑦))
𝐷𝛼2 ( 𝜕 2 )𝛼2 ∕2 𝐾0
−∇
𝜕𝛹 ⎪ 𝐷𝛼1 𝜕2 ⎪ ( 2 )𝛼∕2
𝑖 =⎨ − + − − 𝑛NL (𝐼)⎬ 𝛹 , (5) 1
𝜕𝜁 𝛼1 −1 𝜕𝜉 2 𝛼2 −1 𝜕𝜉 2 𝑛 = 𝑝 + 𝑞2 𝑑𝑝𝑑𝑞
⎪ 𝐾0 𝐾0 0 ⎪ (2𝜋)2 ∬R2 (13)
⎩ ⎭
where 𝜁 is the propagation distance, 𝜉 is the transverse coordinate, 𝐾0 = × 𝜓(𝜉, 𝜂) exp[𝑖𝑝(𝑥 − 𝜉) + 𝑖𝑞(𝑦 − 𝜂)]𝑑𝜉𝑑𝜂.
∬R2
2𝜋∕𝜆 is the propagation constant of the carrier wave with wavelength
Here 𝑝 and 𝑞 are wavenumbers conjugate to coordinates 𝑥 and 𝑦,
𝜆,𝐷𝛼1,2 are relative diffraction coefficients which are determined by
respectively.
the structure of the phase mask(s) which realize the fractional diffrac-
tion [8,10], 𝑛0 is the background refractive index, and 𝐼 = |𝛹 (𝜁 , 𝜉)|2 On the other hand, another natural 2D extension of Eq. (6) for
is the light intensity. The nonlinearity of the setting is represented by anisotropic fractional bulk waveguides, with different LIs acting along
2 the 𝑥 and 𝑦 directions, is defined as
the CQ correction to the √ refractive index, 𝑛NL (𝐼) = 𝑛2 𝐼 + 𝑛4 𝐼 . The [ ( )𝛼1 ∕2 ( )𝛼2 ∕2 ]
substitution 𝜓(𝑧, 𝑥) = ||𝑛4 || ∕ ||𝑛2 ||𝛹0 𝛹 (𝜁 , 𝜉), where 𝛹0 is a character-
−1
𝜕𝜓 𝜕2 𝜕2
√ ( ) ( ) 𝑖 = 𝐷1 − + 𝐷2 − 𝜓 − 𝑔|𝜓|2 𝜓 + 𝛾|𝜓|4 𝜓. (14)
istic amplitude, 𝑥 = 𝐾0 𝑛2 2∕ 𝑛0 ||𝑛4 || 𝜉, and 𝑧 = 𝐾0 𝑛22 ∕ 𝑛0 ||𝑛4 || 𝜁 casts
𝜕𝑧 𝜕𝑥2 𝜕𝑦2
Eq. (5) in a normalized form: The Hamiltonian related to Eq. (14) is
[ ( )𝛼1 ∕2 ( )𝛼2 ∕2 ] { [ ( )𝛼1 ∕2 ( )𝛼1 ∕2 ] }
𝜕𝜓 𝜕2 𝜕2
𝑖 = 𝐷1 − + 𝐷2 − 𝜓 − 𝑔|𝜓|2 𝜓 + 𝛾|𝜓|4 𝜓. (6) 𝐻12 = 𝑑𝑥𝑑𝑦
𝜕2
𝜓 ∗ 𝐷1 − 2
𝜕2
+ 𝐷2 − 2
𝑔 𝛾
𝜓 − |𝜓|4 + |𝜓|6 .
𝜕𝑧 𝜕𝑥2 𝜕𝑥2 ∬R2 𝜕𝑥 𝜕𝑦 2 3

Here the effective positive diffraction coefficients are (15)

√ 2−𝛼1,2
⎛ ⎞
⎜ 𝑛0 ||𝑛4 || ⎟ 2.2. 1D and 2D stationary solutions
𝐷1,2 ≡ 𝐷𝛼1,2 ⎜ (7)
|𝑛2 | ⎟
,
⎜ | | ⎟
⎝ ⎠ The stationary solutions of Eqs. (6), (11), and (14), with real prop-
and the effective nonlinearity coefficients, agation constant 𝑘, are sought for as
( ) ( )
𝑔 ≡ 𝛹02 sgn 𝑛2 , 𝛾 ≡ −𝛹04 sgn 𝑛4 , (8) 𝜓(𝑥, 𝑧) = 𝜙(𝑥; 𝑘)𝑒𝑖𝑘𝑧 , 𝜓(𝑥, 𝑦, 𝑧) = 𝜙(𝑥, 𝑦; 𝑘)𝑒𝑖𝑘𝑧 , (16)

correspond to the cubic self-focusing and quintic defocusing in the cases with confined amplitude functions (lim|𝑥|→∞ 𝜙
of 𝑔 > 0 and 𝛾 > 0, respectively. (𝑥; 𝑘) = 0, lim|𝑥|,|𝑦|→∞ 𝜙(𝑥, 𝑦; 𝑘) = 0) obeying the following equations:
The realization of the fractional diffraction and GVD, proposed and
[ ( )𝛼1 ∕2 ( )𝛼2 ∕2 ]
experimentally created in Refs. [8,10], respectively, suggest that the 𝜕2 𝜕2
scheme may be implemented in the cavity with the longitudinal size −𝑘𝜙 = 𝐷1 − + 𝐷2 − 𝜙 − 𝑔|𝜙|2 𝜙 + 𝛾|𝜙|4 𝜙, (17)
𝜕𝑥2 𝜕𝑥2
∼ 10 cm, using light beams with carrier wavelength 𝜆 ∼ 1 μm and power
∼ 100 kW. The fractional diffraction may be emulated by the hologram in the 1D case corresponding to Eq. (6), or
[ ( )𝛼 ∕2 ( )𝛼 ∕2 ]
with the size ∼ 1 cm, and the expected transverse width of the spatial
−𝑘𝜙 = 𝐷1 −∇2 1 + 𝐷2 −∇2 2 𝜙 − 𝑔|𝜙|2 𝜙 + 𝛾|𝜙|4 𝜙, (18)
solitons will be ∼ 100 μm.
While there are different definitions of fractional derivatives, the [ ( )𝛼1 ∕2 ( )𝛼2 ∕2 ]
𝜕2 𝜕2
one relevant to the realization in fractional quantum mechanics and −𝑘𝜙 = 𝐷1 − + 𝐷2 − 𝜙 − 𝑔|𝜙|2 𝜙 + 𝛾|𝜙|4 𝜙. (19)
optics is the Riesz derivative, based on the juxtaposition of the direct 𝜕𝑥2 𝜕𝑦2
and inverse Fourier transforms, 1 and 1−1 , [1,4,7,53]: in the 2D cases corresponding to Eqs. (11) and (14), respectively.
( )𝛼∕2 Stationary Eqs. (17), (18) and (19) can be numerically solved by
𝜕2 [ ] 1 means of the Newton-conjugate gradient method [54]. Families of the
− 𝜓(𝑥) = 1−1 |𝑝|𝛼 1 (𝜓(𝑥)) ≡ 𝑑𝑝|𝑝|𝛼 𝑑𝜉𝑒𝑖𝑝(𝑥−𝜉) 𝜓(𝜉),
𝜕𝑥2 2𝜋 ∫R ∫R soliton solutions are characterized by their integral power (norm) 𝑁,
(9) defined as

where 𝑝 is the wavenumber conjugate to coordinate 𝑥. 𝑁1D (𝑘) = |𝜙(𝑥; 𝑘)|2 𝑑𝑥, 𝑁2D (𝑘) = |𝜙(𝑥, 𝑦; 𝑘)|2 𝑑𝑥𝑑𝑦, (20)
∫R ∬R2
Note that Eq. (6) is associated with the variational principle 𝑖𝜕𝜓∕
and considered as functions of propagation constant 𝑘.
(𝜕𝑧) = 𝛿𝐻1 ∕(𝛿𝜓 ∗ ), where the Hamiltonian is
[ ( )𝛼1 ∕2 ]
∑ 𝜕2 𝑔 𝛾
𝐻1 = 𝜓 ∗
𝐷𝑗 − 4 6
𝜓 − |𝜓| + |𝜓| 𝑑𝑥, (10) 2.3. Stability analysis of stationary solutions
∫R 𝜕𝑥2 2 3
𝑗=1,2

∗ standing for the complex conjugate. Taking into regard the defi- Once the stationary soliton solutions (16) of Eqs. (6), (11), and (14)
nition of the Riesz derivative (9), it is easy to demonstrate that the are found by numerically solving the stationary Eqs. (17), (18) and
Hamiltonian (10) is real, as it should be. (19), it is necessary to analyze their stability. The linear stability of
1D solitons can be investigated by adding a perturbation to them with
To model the light propagation in the isotropic dual-fractional bulk
an infinitesimal amplitude 𝜖,
medium, Eqs. (6) and (10) can be extended to their 2D versions,
[ ]
𝜕𝜓 [ ( 2 )𝛼1 ∕2 )𝛼 ∕2 ]
∗
( 𝜓(𝑥, 𝑧) = 𝜙(𝑥; 𝑘) + 𝜖(𝑢(𝑥)𝑒𝛿𝑧 + 𝑣∗ (𝑥)𝑒𝛿 𝑧 ) 𝑒𝑖𝑘𝑧 . (21)
𝑖 = 𝐷1 −∇ + 𝐷2 −∇2 2 𝜓 − 𝑔|𝜓|2 𝜓 + 𝛾|𝜓|4 𝜓 (11)
𝜕𝑧
Here 𝑢(𝑥) and 𝑣(𝑥) are perturbation factors, to be found as solutions of
[cf. Eq. (4)] and
the linearized equations,
[ ]
∑ ( )𝛼 ∕2 𝑔 𝛾 ( )( ) ( )
𝐻2 = 𝜓∗ 𝐷𝛼𝑗 −∇2 1 𝜓 − |𝜓|4 + |𝜓|6 𝑑𝑥𝑑𝑦, (12) 𝐿11 𝐿12 𝑢(𝑥) 𝑢(𝑥)
∬R2 2 3 𝑖 ∗ ∗ =𝛿 , (22)
𝑗=1,2 −𝐿12 −𝐿11 𝑣(𝑥) 𝑣(𝑥)

3
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

Fig. 1. (a1) The dependence of propagation constant 𝑘 on power 𝑁for the family of 1D fundamental solitons at 𝑔 = 𝛾 = 1 and 𝛼1 = 2, 𝛼2 = 1.5, as produced by the numerical and
variational solutions of Eq. (17). (a2) The comparison of the numerically found soliton profile (the solid line) and its VA counterpart (the dashed line) at 𝑘 = 0.1. (a3) The stable
evolution of the soliton from panel (a2). (b1) The 1D-soliton’s power 𝑁 vs. quintic coefficient 𝛾 at 𝑔 = 1, 𝑘 = 0.1 and 𝛼1 = 1, 𝛼2 = 1. (b2) The comparison of the numerically found
soliton profile (the solid line) and its VA counterpart (the dashed line) at 𝛾 = 0.04. (b3) The unstable evolution of the soliton from panel (b2). (c1) The 1D-soliton’s power 𝑁 vs.
LI 𝛼2 at 𝑔 = 𝛾 = 1, 𝑘 = 0.1, and 𝛼1 = 2. (c2) The comparison of the numerically found soliton profile (solid line) and its VA counterpart (dashed line) at 𝛼2 = 1. (c3) The stable
evolution of the soliton from panel in (c2). The diffraction coefficients are 𝐷1 = 𝐷2 = 1∕2.

where application is, however, limited. Its accuracy, which depends on the
[ ( )𝛼1 ∕2 ( )𝛼2 ∕2 ] selection of the test function (ansatz), is evaluated below.
𝜕2 𝜕2
𝐿11 ≡ − 𝐷1 − + 𝐷2 − 𝜙 + 2𝑔|𝜙|2 − 3𝛾|𝜙|4 − 𝑘, To apply the VA, we note that Eq. (17) for real 𝜙(𝑥) can be derived
𝜕𝑥2 𝜕𝑥2
from the following Lagrangian:
𝐿12 ≡ 𝑔𝜙2 − 2𝛾|𝜙|2 𝜙2 ,
1
(23) 𝐿= 𝑑𝑝[𝐷1 |𝑝|𝛼1 + 𝐷2 |𝑝|𝛼2 ] 𝑑𝜉𝑑𝑥𝑒𝑖𝑝(𝑥−𝜉) 𝜙(𝑥)𝜙(𝜉)
4𝜋 ∫R ∬R2
[ ] (24)
and 𝛿 ≡ 𝛿𝑅 + 𝑖𝛿𝐼 is an eigenvalue (generally speaking, complex), with 𝑘 𝑔 𝛾
+ 𝑑𝑥 𝜙2 (𝑥) − 𝜙4 (𝑥) + 𝜙6 (𝑥) .
the instability growth rate determined by 𝛿𝑅 , if it is positive. System ∫R 2 4 6
(22) can be solved numerically by means of the Fourier-collocation A simple variational ansatz, approximating soliton solution to Eq. (17),
and Newton-conjugate gradient methods [54]. The results of the linear can be adopted in the form of the Gaussian,
stability analysis should then be verified by direct simulations of the
perturbed propagation. The linear-stability analysis for the 2D Eq. (11) 𝜙(𝑥) = 𝐴 exp[−𝑥2 ∕(2𝑊 2 )], (25)
or (14) can be developed similarly.
with real amplitude
√ 2 𝐴 and width 𝑊 . Power (20) of the 1D ansatz
3. The 1D dual-LI system (25) is 𝑁 = 𝜋𝐴 𝑊 , which may be considered as a variational
parameter, instead of 𝐴. The substitution of ansatz (25) in Lagrangian
In this section, we consider the formation and dynamics of solitons (24) yields the corresponding effective Lagrangian, as a function of the
in the framework of Eq. (6). By means of rescaling, one can fix 𝐷1 = free parameters, 𝑁 and 𝑊 :
𝐷2 = 1∕2 (provided that 𝛼1 ≠ 𝛼2 ) and 𝑔 = +1 in Eq. (6) (assuming that [ ( ) ( ) ]
the cubic term has the focusing sign), keeping the quintic coefficient 𝛾 𝑁 𝛼1 + 1 𝛼2 + 1
𝐿ef f = √ 𝐷1 𝛤 𝑊 −𝛼1 + 𝐷2 𝛤 𝑊 −𝛼2
as a freely varying parameter. 2 𝜋 2 2
(26)
𝑘 𝑔 𝛾
+ 𝑁 − √ 𝑁 2 𝑊 −1 + √ 𝑁 3 𝑊 −2 ,
3.1. The variational approximation (VA) 2 4 2𝜋 6 3𝜋

Localized solutions of Eq. (17) can be sought for in an approximate where 𝛤 (⋅) is the Gamma-function. Then, for given 𝑘, the values of 𝑁
analytical form by means of VA [13,50,51]. It is a useful tool, whose and 𝑊 can be obtained from the Euler–Lagrange equations, 𝜕𝐿ef f ∕𝜕𝑁 =

4
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

𝜕𝐿ef f ∕𝜕𝑊 = 0, viz., TSs in the 2D cubic NLSE with the normal diffraction (𝛼 = 2) [34,35].
( ) ( ) A characteristic feature of the TS family in any model is its degeneracy,
𝛼1 + 1 𝛼2 + 1 √ 𝑔 𝑁
𝐷1 𝛤 𝑊 −𝛼1 + 𝐷2 𝛤 𝑊 −𝛼2 + 𝑘 𝜋 − √ in the sense that all the solitons, with any value of the propagation
2 2 2𝑊 constant, have a single (critical) value of the power, 𝑁crit . For the sake
𝛾 𝑁 2
+√ = 0, (27a) of the completeness, we briefly reproduce here the derivation of the
2
3𝜋 𝑊 condition for the possibility of the critical or super-critical collapse
( ) ( )
𝛼1 + 1 𝛼2 + 1 √ in the framework of Eq. (6). To this end, note that the respective
𝐷1 (𝛼1 − 1)𝛤 𝑊 −𝛼1 + 𝐷2 (𝛼2 − 1)𝛤 𝑊 −𝛼2 − 𝑘 𝜋 Hamiltonian with 𝐷1 = 𝐷2 = 1∕2, 𝛼1 = 𝛼2 = 1 is
2 2
[ ( )1∕2 ]
𝑔 𝑁 𝛾 𝑁2
+ √ − √ = 0. (27b) 𝜕2 𝑔 𝛾
𝑊 2 𝐻1 = 𝜓∗ − 𝜓 − |𝜓|4 + |𝜓|6 𝑑𝑥
2 2 3 3𝜋 𝑊 ∫R 𝜕𝑥2 2 3 (31)
Eqs. (27a) and (27b) can be solved numerically by means of the ≡ 𝐻diffraction + 𝐻focusing + 𝐻defocusing .
relaxation method.
Considering a localized state with radius 𝐿 and amplitude 𝐵 [e.g., 𝜙(𝑥)
( )
3.2. 1D solitons = 𝐵 exp −𝑥2 ∕(2𝐿2 ) , as in Eq. (25)], an obvious estimate for the power
is 𝑁 ∼ 𝐵 2 𝐿. Similarity, the scaling of the diffraction, self-focusing cubic
First, we address a relation between the soliton’s power 𝑁 and and self-defocusing quintic terms in the Hamiltonian (31), following the
propagation constant 𝑘, as produced by the VA and numerical solutions. variation of 𝐿, can be estimated as:
A typical example, for a set of parameter 𝛼1 = 2, 𝛼2 = 1.5, 𝑔 = 𝛾 = 1, is
exhibited in Fig. 1(a1), which demonstrates perfect agreement between 𝐻diffraction ∼ 𝑁𝐿−1 , 𝐻focusing ∼ −𝑔𝑁 2 𝐿−1 , 𝐻defocusing ∼ 𝛾𝑁 3 𝐿−2 .
the VA and numerical results. The 𝑘(𝑁) curve displayed in this figure (32)
complies with the Vakhitov–Kolokolov (VK) criterion, 𝑑𝑘∕𝑑𝑁 > 0,
which is a necessary condition for stability of solitons supported by a For a fixed power 𝑁, in the absence of the self-defocusing quintic
dominant self-focusing nonlinearity [33,55]. The linear-stability results term (𝛾 = 0), the collapse occurs if |𝐻focusing | diverges faster than
and systematic numerical simulations demonstrate that all the solitons 𝐻diffraction at 𝐿 → 0. Thus, the critical or super-critical collapse, which
are stable at 𝛾 > 0.04 (the instability of the solitons at very small 𝛾 is sets in, severally, at 𝑁 > 𝑁crit , or for arbitrarily small 𝑁, takes place
discussed below). An illustration of the dynamical stability is presented at 𝛼1 = 𝛼2 = 1 or 𝛼1 = 𝛼2 < 1, respectively.
in Fig. 1(a3).
In Fig. 1(a1), the cutoff value of propagation constant 𝑘 for the 3.3.1. The self-defocusing quintic nonlinearity
existence of the solitons is 𝑘max = 3∕16, which is a known feature The self-defocusing quintic term (𝛾 > 0) arrests the onset of the
of the CQ nonlinearity, irrespective of the spatial dimension [22,37– collapse. The dependence of power 𝑁 on coefficient 𝛾 for the 1D-soliton
40,42,56–60]. In particular, Eq. (6) with 𝛼1 = 𝛼2 = 2 and 𝑔 = 𝛾 = 1 family with 𝛼1 = 𝛼2 = 1, 𝑔 = 1 and 𝑘 = 0.1 is displayed in Fig. 1(b1).
gives rise to the well-known family of 1D bright solitons [61], The evolution of the soliton profile with the increase of 𝛾 is displayed in
√
( ) √ Fig. 2(b). The growth of the amplitude and width of the stable solitons
𝜓sol CQ = 2 exp(𝑖𝑘𝑧)√
𝑘
√ √ ( √ ), (28) with 𝛾, observed in this figure, may be qualitatively explained by the
1 + 1 − 16𝑘∕3 cosh 2 𝑘𝑥
comparison with the exact soliton solution for 𝛼1 = 𝛼2 = 2, given by
in which the propagation constant 𝑘 takes values in the interval of Eq. (28) (in which 𝛾 was absorbed by rescaling of the equation).
0 < 𝑘 < 𝑘max ≡ 3∕16. In the limit case 𝑘 → 3∕16, bright soliton (28) As mentioned above, it is found that the soliton is unstable only
goes over into a stationary front solution [62] at very small values of 𝛾 ≤ 0.04, in agreement with the expectation
√ that the dominant cubic self-focusing produces the destabilizing effect
( ) ( ) 3∕4
3𝑖 (at so small values of 𝛾, the results are affected by a finite size of the
𝜓front CQ = exp 𝑧 √ , (29)
16 1 + exp(± 3𝑥) integration domain and finite propagation distance). The profile and
evolution of the soliton at 𝛾 = 0.04, which still seems unstable, are
which connects asymptotic flat states, 𝜓 = 0 and
√ exhibited in Figs. 1(b2) and (b3), respectively. The instability does not
𝜓 = ( 3∕2) exp (3𝑖𝑧∕16) . (30) destroy the soliton, but transforms it into a breather.
The impact of LI 𝛼2 for fixed 𝛼1 = 2, 𝑔 = 𝛾 = 1 and 𝑘 = 0.1 is
Accordingly, the solitons develop a flat-top shape in the limit of 𝑘 →
explored too. Power 𝑁 versus the LI 𝛼2 is shown in Fig. 1(c1), which
𝑘max , and the relaxation method, if applied for solving system (27),
demonstrates good accuracy of the VA. The linear-stability analysis
diverges at 𝑘 > 0.1765, as the Gaussian ansatz cannot adequately
and direct propagation corroborate that all solitons are stable when 𝛼2
approximate the flat-top solitons at 𝑘 → 𝑘max . An example of a soliton
varies between 1 to 2. It is also seen in Fig. 2(c) that the amplitude
with 𝑘 = 0.1 is exhibited in Fig. 1(a2), showing that the VA prediction
and width of the 1D solitons decrease with the increase of 𝛼2 . A typical
matches well with the numerical solution, except for inability of the
profile and stable evolution of the soliton for 𝛼2 = 1 are exhibited in
Gaussian ansatz to reproduce the power-law decaying tails of the
Figs. 1(c2) and (c3), respectively.
soliton [53,63]. The evolution of the numerically found shape of the
1D solitons, following the increase of 𝑘, is displayed in Fig. 2(a). It
demonstrates the transition to the flat-top shape of the solitons at 3.3.2. The self-focusing quintic nonlinearity
𝑘 → 𝑘max . Next, we address the case of 𝑔 = 1, 𝛾 = −1, where both the cubic
and quintic terms are self-focusing. It is known that, in the case of
3.3. The effect of the quintic coefficient 𝛾 on the 1D solitons the normal diffraction (𝛼1 = 𝛼2 = 2), while all one-dimensional TSs,
obtained for 𝑔 = 0 and 𝛾 < 0, were unstable [42,49], in the case of
Next, we consider the effect of variation of coefficient 𝛾 in front 𝑔 > 0 and 𝛾 < 0 the exact soliton solutions of the 1D CQ-NLS equation
of the quintic term in Eq. (6), while the strength of the self-focusing are stable against small perturbations. However, once the Lévy index
cubic term is fixed a set above, 𝑔 = +1. Recent results [50,51,64] reveal in the “single-LI’’ FNLSE is set to be 𝛼 < 2, the solitons are unstable
that the FNLSE (6) gives rise to quasi-TSs for 𝛼1 = 𝛼2 = 1 and 𝛾 = 0, when propagation constant 𝑘 is greater than some threshold value.
which are made unstable by the possibility of thecritical collapse in the Hence, following the ideas mentioned above, we consider the option
same case [32,33], similar to the situation for the commonly known to include an additional diffraction term into the “single-LI’’ FNLSE,

5
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

Fig. 2. The evolution of profiles of the 1D solitons, following the variation of propagation constant 𝑘 (a), quintic coefficient 𝛾 (b), or LI 𝛼2 (c), is displayed by means of heat
maps of 𝜙(𝑥), see Eq. (16). In panel (a), the range of 0 < 𝑘 < 𝑘max is shown, outside of which solitons do not exist. The parameters are the same as Fig. 1.

Fig. 3. (a1) The dependence of the 1D fundamental-soliton’s power 𝑁 on propagation constant 𝑘 at 𝑔 = 1, 𝛾 = −1, 𝐷2 = 1∕2, and 𝛼1 = 2, 𝛼2 = 1.5, as produced by the numerical
solutions of Eq. (17), which corresponds to the case when both the cubic and quintic terms are self-focusing. The red and purple lines correspond, severally, to 𝐷1 = 1∕2 or 𝐷1 = 0.
The solid and dashed lines represent stable and unstable solitons, respectively. Labels (a2, a3) correspond to the solitons displayed in panels (a2, a3). (a2) The unstable evolution
of the soliton for 𝑘 = 1 and 𝐷1 = 0. (a3) The stable evolution of the soliton for 𝑘 = 1 at 𝐷1 = 1∕2, in the situation when the interplay of the fractional diffraction with 𝛼2 = 1.5
and quintic self-focusing would lead, by itself, to the instability driven by the super-critical collapse. The insets display the corresponding initial soliton profiles.

and then check the effect of this on the soliton stability. To this end, 4.1. 2D fundamental isotropic solitons
we first set
Proceeding to the 2D setting, we first address the isotropic system
𝐷1 = 0, 𝐷2 = 1∕2, 𝛼2 = 1.5
based on Eq. (11). The respective equation (18) can be derived from
in Eq. (6), which amounts to the “single-LI’’ FNLSE with 𝛼 = 1.5. the following Lagrangian for real 𝜙(𝑥, 𝑦):
As indicated in Fig. 3(a1), in this case the solitons are stable only at 𝐷 [( )𝛼 ∕2 ]
𝑘 ≤ 0.4, there being two soliton branches with opposite signs of the 𝐿= 1 −∇2 1 𝜙(𝑥, 𝑦) 𝜙(𝑥, 𝑦)𝑑𝑥𝑑𝑦
2 ∬R2
slopes, 𝑑𝑁∕𝑑𝑘, the respective (in)stability being consistent with the 𝐷2 [( )𝛼 ∕2 ]
𝑘
VK criteria. A typical soliton profile and the corresponding unstable + −∇2 2 𝜙(𝑥, 𝑦) 𝜙(𝑥, 𝑦)𝑑𝑥𝑑𝑦 + 𝜙2 (𝑥, 𝑦)𝑑𝑥𝑑𝑦
2 ∬R2 2 ∬R2
evolution are displayed in Fig. 3(a2). Next, we include the additional 𝑔 𝛾
diffraction term into the FNLSE, setting − 𝜙4 (𝑥, 𝑦)𝑑𝑥𝑑𝑦 + 𝜙6 (𝑥, 𝑦)𝑑𝑥𝑑𝑦, (35)
4 ∬R2 6 ∬R2
𝐷1 = 𝐷2 = 1∕2, 𝛼1 = 2, 𝛼2 = 1.5 (33) ( )𝛼 ∕2
where −∇2 1,2 𝜙(𝑥, 𝑦) is defined as per Eq. (13). Similar to the 1D
in Eq. (6), i.e., adding the non-fractional diffraction term (with 𝛼1 = 2). ansatz (25), the Gaussian may be adopted as the variational ansatz for
Systematic numerical simulations demonstrate that the soliton-stability 2D solitons:
( 2 )
region expands, as expected, see Fig. 3(a1). We also display the soliton 𝑥 + 𝑦2
𝜙(𝑥, 𝑦) = 𝐴 exp − , (36)
profile with 𝑘 = 1, and its stable evolution, in Fig. 3(a3). Thus, 2𝑊 2
an essential conclusion is that the synergetic effect of two different
with power 𝑁 = 𝜋𝐴2 𝑊 2 . The substitution of this ansatz in La-
diffraction terms may produce stable 1D solitons in the situation when
grangian (35) yields the corresponding effective Lagrangian, which can
the interplay of the single fractional-diffraction term (with LI 𝛼2 = 1.5)
be conveniently written as a function of 𝑊 and 𝑁, cf. Eq. (26):
and quintic self-focusing would produce the strong instability driven by ( ) ( )
the super-critical collapse. 𝐷 𝛼1 + 2 𝐷 𝛼2 + 2
𝐿ef f = 1 𝛤 𝑁𝑊 −𝛼1 + 2 𝛤 𝑁𝑊 −𝛼2
2 2 2 2
(37)
4. The 2D dual-LI CQ-FNLS systems 𝑘 𝑔 2 −2 𝛾
+ 𝑁− 𝑁 𝑊 + 𝑁 3 𝑊 −4 .
2 8𝜋 18𝜋 2
In this section, we consider the formation and stability of fundamen- Then, for given 𝑁, the value of width 𝑊 is predicted by the Euler–
tal and vortex solitons, with the winding numbers 𝑆 = 0 and 1, 2, 3, in Lagrange equation, 𝜕𝐿ef f ∕𝜕𝑊 = 0:
the framework of 2D Eqs. (11) and (14). Once again making use of the ( ) ( )
scaling freedom, we set 𝐷1 = 𝐷2 = 1∕2 in Eq. (18) and 𝐷1 = 𝐷2 = 1 in 𝐷 𝛼1 + 2 𝐷 𝛼2 + 2 𝑔
𝛼1 1 𝛤 𝑊 −𝛼1 −1 + 𝛼2 2 𝛤 𝑊 −𝛼2 −1 − 𝑁𝑊 −3
Eq. (19). We also fix 2 2 2 2 4𝜋
2𝛾 2 −5
𝛼1 = 2, 𝛼2 = 1.5, 𝑔=𝛾=1 (34) + 𝑁 𝑊 = 0.
9𝜋 2
in most cases, except for the vortex solitons produced by Eq. (19). (38)

6
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

Fig. 4. (a) Power 𝑁 of the 2D isotropic fundamental soliton vs. the propagation constant, 𝑘, as produced by the numerical solution of Eq. (18) and the respective VA, for
𝐷1 = 𝐷2 = 1∕2 and other parameters fixed as per Eq. (34). The red dashed and solid lines and the purple dotted one represent, respectively, unstable and stable numerically
found solitons, and their VA-predicted counterparts. (b) Cross-section profiles 𝜙 (𝑥, 0) of the numerically found 2D solitons corresponding to the propagation constant taking values
0.001 ≤ 𝑘 ≤ 0.1855. (c) and (e): Profiles 𝜙(𝑥, 0) for 𝑘 = 0.0145 and 𝑘 = 0.0435, respectively, as obtained from the numerical and VA solutions. (d) and (f): The unstable and stable
propagation of the (weakly perturbed) solitons from panels (c) and (e), respectively, displayed by means of isosurfaces of the local power |𝜙|2 .

The relation between power 𝑁 and propagation constant 𝑘 for the 4.2. Fundamental solitons in the system with anisotropic fractional diffrac-
2D solitons, produced by the numerical solution and its VA counter- tion
part derived from Eq. (38) and the second Euler–Lagrange equation,
𝜕𝐿ef f ∕𝜕𝑁 = 0, is presented in Fig. 4(a). The VA procedure diverges at Next, we consider fundamental solitons supported by the fractional-
𝑘 > 0.1175, similar to what is mentioned above for the 1D case. It is seen diffraction terms acting along the 𝑥 and 𝑦 direction with different LIs,
that the 𝑘(𝑁) dependence is composed of top and bottom branches, as defined by Eqs. (14) and (19). One may then expect that the model
which exist at powers exceeding a minimum (threshold) value, gives rise to elliptically shaped solitons, with the ellipse’s major axis
aligned with the direction corresponding to smaller LI.
𝑁thr ≈ 18.5018, (39) The Lagrangian of Eq. (19) for real 𝜙(𝑥, 𝑦) is (cf. expression (35))
𝐷 [( )𝛼 ∕2 ]
with the positive and negative slopes, 𝑑𝑘∕𝑑𝑁 > 0 and 𝑑𝑘∕𝑑𝑁 < 0, 𝐿= 1 −𝜕𝑥𝑥 1 𝜙(𝑥, 𝑦) 𝜙(𝑥, 𝑦)𝑑𝑥𝑑𝑦
respectively. In accordance with the above-mentioned VK criterion [33, 2 ∬R2
𝐷2 [( )𝛼 ∕2 ]
𝑘
55], the linear-stability analysis and systematic simulations of the + −𝜕𝑦𝑦 2 𝜙(𝑥, 𝑦) 𝜙(𝑥, 𝑦)𝑑𝑥𝑑𝑦 + 𝜙2 (𝑥, 𝑦)𝑑𝑥𝑑𝑦
perturbed evolution of the 2D solitons demonstrate that the top and 2 ∬R2 2 ∬R2
bottom branches are, respectively, stable and unstable ones. Note that, 𝑔 𝛾
− 𝜙4 (𝑥, 𝑦)𝑑𝑥𝑑𝑦 + 𝜙6 (𝑥, 𝑦)𝑑𝑥𝑑𝑦, (41)
in the case of the normal diffraction (𝛼1 = 𝛼2 = 2), the existence of 4 ∬R2 6 ∬R2
two branches with the opposite stability is a characteristic feature of ( )𝛼 ∕2 ( )𝛼 ∕2
where −𝜕𝑥𝑥 1 𝜙(𝑥, 𝑦) and −𝜕𝑦𝑦 2 𝜙(𝑥, 𝑦) are defined as per Eq. (9).
the 3D version of the NLSE with the CQ nonlinearity [40]. A similar The respective variational ansatz is adopted in the form of an
feature was recently reported for the single-LI 2D FNLSE, namely, the anisotropic Gaussian,
one with 𝛼1 = 𝛼2 < 2, and the CQ nonlinearity [60]. The similarity is ( )
explained by the fact the 2D FNLSE with LI < 2 and cubic self-focusing 𝑥2 𝑦2
𝜙(𝑥, 𝑦) = 𝐴 exp − − , (42)
gives rise to the supercritical collapse, as well as its 3D counterpart with 2𝑊12 2𝑊22
the normal diffraction (LI = 2) [42].
with unequal widths 𝑊1 and 𝑊2 along the 𝑥 and 𝑦 directions, and the
In terms of the VA, 𝑁thr may be identified as a point at which respective norm 𝑁 = 𝜋𝐴2 𝑊1 𝑊2 , unlike the above ansatz (36). The
two roots of Eq. (38) for 𝑊 merge into one, which is determined by substitution of expression (42) in Lagrangian given by Eq. (41) yields
equation 𝜕 2 𝐿ef f ∕𝜕𝑊 2 = 0 [42]. In the combination with Eq. (38), it the corresponding effective Lagrangian:
yields ( ) ( )
𝐷 𝛼1 + 1 −𝛼 𝐷 𝛼2 + 1 −𝛼
𝐿ef f = √1 𝛤 𝑁𝑊1 1 + √2 𝛤 𝑁𝑊2 2
(VA)
𝑁thr ≈ 19.2050, (40) 2 𝜋 2 2 𝜋 2
𝑘 𝑔 2 −1 −1 𝛾
+ 𝑁− 𝑁 𝑊1 𝑊2 + 𝑁 3 𝑊1−2 𝑊2−2 . (43)
which is close to its numerical counterpart given by Eq. (39). 2 8𝜋 18𝜋 2
Profiles of the cross-section 𝜙(𝑥, 0) for the 2D solitons with propaga- Then, for given 𝑘, values of 𝑁, 𝑊1 and 𝑊2 are obtained from the
tion constant taking values 0.001 ≤ 𝑘 ≤ 0.1855 are displayed in Fig. 4(b), corresponding Euler–Lagrange equations, 𝜕𝐿ef f ∕𝜕𝑁 = 𝜕𝐿ef f ∕𝜕𝑊1 =
which features the flat-top shape at 𝑘 → 𝑘max with 𝜙(𝑥, 𝑦) approaching 𝜕𝐿ef f ∕𝜕𝑊2 = 0. The solution of these equations diverge at 𝑘 > 0.12,
the above-mentioned
√ largest value admitted by 1D solution (28), viz., when ansatz (42) and Lagrangian (43) are irrelevant, similar to the
𝜙max = 3∕2, see Eq. (30). The solitons pertaining to 𝑘 = 0.0145 and above-mentioned divergence in the isotropic 2D model.
0.0435, which are, respectively, unstable and stable ones, are presented Results for the fundamental solitons produced by the numerical and
in Figs. 4(c) and (e). The corresponding unstable and stable evolution variational solutions of Eq. (19) with 𝑔 = 𝛾 = 1 and 𝛼1 = 2, 𝛼1 = 1.5
is exhibited, severally, in Figs. 4(d) and (f). are summarized in Fig. 5(a1), in the form of dependence 𝑘(𝑁). It is

7
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

Fig. 5. (a1) Power 𝑁 of the 2D anisotropic fundamental soliton vs. the propagation constant, 𝑘, as obtained from the numerical and variational solutions of Eq. (19), for 𝐷1 = 𝐷2 = 1
and other parameters fixed as per Eq. (34). (a2) and (a3): The evolution of the cross-section profiles of the anisotropic solitons, 𝜙(𝑥, 0) and 𝜙(0, 𝑦) for the propagation constant
taking values 0.001 ≤ 𝑘 ≤ 0.1823. (b1) and (c1): Unstable and stable evolution of the (weakly perturbed) 2D anisotropic solitons corresponding to 𝑘 = 0.0160 and 0.1800, respectively.
(b2), (b3) and (c2), (c3): Comparison of the initial and final shapes of the evolving solitons for the same cases.

Fig. 6. The critical value produced by numerical and variational solutions of Eq. (19) at (a) 𝛼2 = 1, (b) 𝛼2 = 1.5 and (c) 𝛼2 = 2, for 𝛼1 varying from 1 to 2. Other parameters are
same as Fig. 5.

seen that, similar to what is displayed in Fig. 4(a), two branches with branches in Fig. 5(a1) are stable and unstable, respectively, in full ac-
positive and negative slopes merge at the critical point, with 𝑘 = 0.0235, cordance with the VK criterion. In direct simulations, unstable solitons
in Fig. 5(a1). In Figs. 5(a2) and (a3), cross-sections 𝜙(𝑥, 0) and 𝜙(0, 𝑦) spontaneously evolve into anisotropic breathers, as shown in Fig. 5(b1).
show the flat-top shapes of the solitons at 𝑘 → 𝑘max . As conjectured The corresponding unstable and stable evolution, and final shapes of
above, the elliptical solitons indeed exhibit their major axis aligned the evolving solitons, are displayed, severally, in Figs. 5(b1, c1) and
with the direction corresponding to the smaller LI. Typical profiles Figs. 5(b3, c3).
of the anisotropic solitons with 𝑘 = 0.0160 and 𝑘 = 0.1800, which Note that the critical (threshold) values of the norm are crucially
belong, respectively, to the unstable and stable branches, are presented important for the existence and stability of fundamental solitons in
in Fig. 5(b2) and (c2). the anisotropic fractional model based on Eq. (19). Therefore, the
As above, the linear stability analysis and direct simulations demon- dependence of numerical and VA-predicted values of 𝑁thr on 𝛼1 with
strate that the 2D-soliton families corresponding to the top and bottom fixed 𝛼2 are presented in Fig. 6, which demonstrates good agreement

8
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

Fig. 7. (a1–a3) Dependences 𝑘(𝑁) for 2D vortex solitons with winding numbers 𝑆 = 1, 2, 3, as produced by the numerical and variational solutions of Eq. (18) for 𝐷1 = 𝐷2 = 1∕2 and
other parameters fixed as per Eq. (34) (all results reported in this section were produced for Eq. (11) with these values of the parameters). (b1–b3) The variation of cross-section
profiles |𝜙(𝑥, 0)| of the vortex solitons with 𝑆 = 1, 2, 3, following the increase of the propagation constant 𝑘.

between the numerical and variational results. It is seen that 𝑁thr with the real-valued amplitude 𝐴 and width 𝑊 . The norm of ansatz
decreases as 𝛼1 increases (for fixed 𝛼2 ). This trend may be explained (45) is
by the fact that the increase of 𝛼1 makes the diffraction stronger in the
system, thus requiring stronger self-focusing nonlinearity to balance it. 𝑁𝑆 = 𝜋𝑆!𝐴2 𝑊 2(1+𝑆) ,
The latter condition is provided by the attenuation of the defocusing √
√ |𝜙(𝑥, 𝑦)|max = 𝐴 𝑆𝑊 𝑒
and√its maximum value, −𝑆∕2 , is attained at
effect of the quintic term with the decrease of 𝑁. A similar conclusion
𝑟 = 𝑆𝑊 , where 𝑟 = 𝑥2 + 𝑦2 .
drawn from the present findings is that, for fixed 𝛼1 , 𝑁thr decreases
The substitution of ansatz (45) in Lagrangian (44) yields the corre-
with the increase of 𝛼2 .
sponding effective Lagrangian, which is again written as a function of
𝑊 and 𝑁, cf. Eq. (37):

4.3. Vortex solitons ( ) ( )

𝛤 𝑆 + 1 + 𝛼1 ∕2 −𝛼1
𝛤 𝑆 + 1 + 𝛼2 ∕2
𝐿ef f = 𝜋𝐷1 𝑁𝑆 𝑊 + 𝜋𝐷2 𝑁𝑆 𝑊 −𝛼2
𝑆! 𝑆!
The possibility to construct 2D vortex solitons in models based 𝜋𝑔 (2𝑆)! 2 −2 𝜋𝛾 (3𝑆)! 3 −4
+𝑘𝑁𝑆 − 𝑁 𝑊 + 𝑁 𝑊 . (46)
on 2D FNLSEs is a subject of obvious interest [13,14,21]. We here 22𝑆+2 (𝑆!)2 𝑆 33𝑆+2 (𝑆!)3 𝑆
consider vortex solitons in the isotropic system with two different For given 𝑘, values of 𝑁𝑆 and 𝑊 are produced by the Euler–Lagrange
fractional-diffraction terms, and in the anisotropic system with the equations, 𝜕𝐿ef f ∕𝜕𝑁𝑆 = 𝜕𝐿ef f ∕𝜕𝑊 = 0. Note that expression (46)
diffraction operators with different LIs, acting along two transverse amounts to Eq. (37) for 𝑆 = 0 (up to a scale factor 1∕2).
coordinates. The corresponding models are based on Eqs. (11) and (14),
Dependencies of the propagation constant 𝑘 on power 𝑁 for the vor-
respectively.
tex solitons with 𝑆 = 1, 2, and 3 are presented, severally, in Figs. 7(a1,
Unlike the real stationary field 𝜙 (𝑥, 𝑦) for the fundamental solitons, a2, a3), in which the stability and instability of the numerically found
solutions for the vortex ones are represented by complex functions. families is designated by the red solid and dashed lines, respectively.
Accordingly, Lagrangian (35) for the 2D real field is replaced by Similar to the fundamental-soliton families, vortex solitons with all
the following Lagrangian corresponding to Eq. (19) for the complex values of 𝑆 populate the top and bottom branches of the 𝑘(𝑁) curves.
function 𝜙(𝑥, 𝑦): The critical powers at which the two branches merge, and below which
[( )𝛼 ∕2 ] no vortex solitons exist, are (as found from the numerical solution)
𝐿 = 𝐷1 −∇2 1 𝜙(𝑥, 𝑦) 𝜙∗ (𝑥, 𝑦)𝑑𝑥𝑑𝑦
∬R2 (𝑆=1) (𝑆=2) (𝑆=3)
[( ] 𝑁thr = 77.12, 𝑁thr = 143.99, 𝑁thr = 212.71. (47)
)𝛼 ∕2
+𝐷2 −∇2 2 𝜙(𝑥, 𝑦) 𝜙∗ (𝑥, 𝑦)𝑑𝑥𝑑𝑦 + 𝑘 |𝜙(𝑥, 𝑦)|2 𝑑𝑥𝑑𝑦
∬R2 ∬R2 Note that these values demonstrate a nearly linear dependence on 𝑁,
𝑔 𝛾 similar to an approximate analytical result obtained in Ref. [66] for
− |𝜙(𝑥, 𝑦)|4 𝑑𝑥𝑑𝑦 + |𝜙(𝑥, 𝑦)|6 𝑑𝑥𝑑𝑦, (44)
2 ∬R2 3 ∬R2 the critical norms of 2D Townes solitons with embedded vorticity (in
( )𝛼 ∕2 the case of the normal diffraction, 𝛼 = 2). It is also relevant to note
with −∇2 𝑗 𝜙(𝑥, 𝑦) defined as per Eq. (13). A straightforward vari- that values (47), obtained for 𝛼1 = 2 and 𝛼2 = 1.5 [see Eq. (34)], are
ational ansatz for vortex solitons with topological charge (winding naturally located between those for 𝛼1 = 𝛼2 = 2 and 𝛼1 = 𝛼2 = 1.5, cf.
number) 𝑆 = 1, 2, 3, …, if written in the Cartesian coordinates, is (cf. Ref. [22].
Ref. [65]) Unlike the situation for 𝑆 = 0, the VK criterion is not sufficient for
( 2 ) the stability of the vortices. Accordingly, numerical results demonstrate
𝑥 + 𝑦2
𝜙(𝑥, 𝑦) = 𝐴(𝑥 + 𝑖𝑦)𝑆 exp − , (45) that only relatively small parts of the upper branches are stable in
2𝑊 2

9
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

Fig. 8. The stable evolution of perturbed vortex solitons produced by Eq. (18), with different winding numbers: 𝑆 = 1, 𝑘 = 0.1500 [panels (a1–a3)]; 𝑆 = 2, 𝑘 = 0.1705 [panels
(b1–b3)]; 𝑆 = 3, 𝑘 = 0.1765 [panels (c1–c3)]. The stable propagation is plotted in the first column. The second and third columns display, severally, the intensity and phase patterns
of the initial states.

Figs. 7(a1, a2, a3), in intervals In particular, the vortex soliton with 𝑆 = 1 and 𝑘 = 0.1380, which is
located very close to the stability boundary, with a small instability
⎧ 0.1390 ≤ 𝑘 < 3∕16 for 𝑆 = 1,
⎪ growth rate, viz., max(𝛿𝑅 ) = 0.0065, splits in two fragments, as shown
⎨ 0.1645 ≤ 𝑘 < 3∕16 for 𝑆 = 2, (48) in Fig. 9(a1). The corresponding initial and final intensity patterns and
⎪ 0.1765 ≤ 𝑘 < 3∕16 for 𝑆 = 3. the initial phase pattern are presented in Figs. 9(a2, a4) and (a3),
⎩
respectively.
Cross-section profiles |𝜙(𝑥, 0)| of the ring-shaped vortex solitons For the vortex soliton with 𝑆 = 2, 𝑘 = 0.0030, and max(𝛿𝑅 ) = 0.1019,
with 𝑆 = 1, 2, 3 for varying propagation constant 𝑘 are displayed in which splits relatively fast in four fragments [see Figs. 9 (b1) and (b4)],
Fig. 7(b1–b3). They also exhibit the trend for the formation of the flat- the corresponding initial intensity and phase patterns are displayed in
top profile, with a relatively small vorticity-carrying hole at the pivot, Figs. 9(b2) and (b3). Finally, an example of the complex development
in the limit of 𝑘 → 𝑘max = 3∕16. of the instability of a vortex with 𝑆 = 3, 𝑘 = 0.0600, and max(𝛿𝑅 ) =
Examples of stable dynamics of the vortex solitons with different 0.1061, which leads to the splitting in a set of six fragments, is displayed
values of the winding number 𝑆 are collected in Fig. 8. For 𝑆 = 1, the in Figs. 9(c1) and (c4). The corresponding initial intensity and phase
intensity and phase pattern for the soliton with 𝑘 = 0.1500 are displayed patterns are plotted in Figs. 9(c2) and (c3).
in Figs. 8(a2) and (a3), while its stable evolution under the action of The number of fragments into which vortex solitons split under the
random-noise perturbations at the 2% amplitude level is presented in action of the azimuthal instability, is often, but not always, equal to
Fig. 8(a1). Fig. 8(b1) demonstrates stable evolution of a vortex soliton integer azimuthal index 𝐽 of the instability eigenmode with the largest
with 𝑆 = 2 and 𝑘 = 0.1705, whose amplitude and phase patterns growth rate [37]. To elucidate the structure of the azimuthal instability
are displayed in Figs. 8(b2) and (b3). Finally, the stable evolution of in this case, the perturbed solutions, with an infinitesimal amplitude 𝜖,
a broad (flat-top) soliton with 𝑆 = 3 at 𝑘 = 0.1765 is presented in should be written in cylindrical coordinates (𝑟, 𝜃, 𝑧) as
Fig. 8(c1), while its amplitude and phase patterns are displayed in [ ∗
]
𝜓(𝑟, 𝜃, 𝑧) = 𝑒𝑖(𝑘𝑧+𝑆𝜃) 𝜙(𝑟) + 𝜖𝑝(𝑟)𝑒+𝑖𝐽 𝜃 𝑒𝛿𝑧 + 𝜖𝑞 ∗ (𝑟)𝑒−𝑖𝐽 𝜃 𝑒𝛿 𝑧 , (49)
Figs. 8(c2) and (c3).
On the other hand, examples of the evolution of unstable vortex where 𝑝(𝑟) and 𝑞(𝑟) are form-factors of the perturbation eigenmode,
solitons are reported in Fig. 9, showing the symmetry-breaking splitting cf. Eq. (21). The respective linearized equations [cf. Eq. (22) in the
instability, which is a generic feature of unstable vortices [37–39,67]. 1D model] were solved in the Cartesian coordinates by means of

10
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

Fig. 9. The unstable evolution (splitting) of vortex solitons: 𝑆 = 1, 𝑘 = 0.1380 [panels (a1–a4)]; 𝑆 = 2, 𝑘 = 0.0300 [panels (b1–b4)]; 𝑆 = 3, 𝑘 = 0.0600 [panels (c1–c4)]. The first
column displays the unstable propagation. The second and third columns depict the intensity and phase patterns of the initial states. The final intensity patterns are plotted in the
fourth column.

Fig. 10. Perturbation eigenmodes for the unstable vortex solitons from Fig. 9. The first, second and third columns relate to the vortex solitons exhibited in the first, second and
third rows of Fig. 9. Intensity and phase patterns of the eigenmodes are plotted in panels (a1–a3) and (b1–b3), respectively.

the Newton-conjugate-gradient method [54]. The eigenmodes with Vortex solitons have also been found in the anisotropic model
the largest growth rate, corresponding to 𝐽 = 2, 4, 6 (see Fig. 10), based on Eq. (19). Naturally, anisotropic (azimuthally deformed) vortex
which pertain to unstable vortex solitons presented in Fig. 9, indeed solitons are more vulnerable to splitting. For this reason, we have found
dominate the development of the splitting instability, as corroborated only a small stability interval for the vortex solitons with topological
by a detailed numerical analysis. charge 𝑆 = 1, see Fig. 11(a1). In particular, the dependence of the

11
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

Fig. 11. (a1) Dependences 𝑁(𝛼2 ) for the 2D vortex solitons with winding numbers 𝑆 = 1, as produced by the numerical solutions of Eq. (19) for 𝐷1 = 𝐷2 = 1 and other parameters
fixed as per Eq. (50). Label (a2) corresponds to the soliton presented in panel (a2). (a2–a3) The intensity and phase patterns of the initial state marked in (a1) with 𝛼2 = 1.94.
(a4–a5) The intensity of the vortex soliton (a2) at 𝑧 = 100 and 700 in the course of the propagation. (a6) The stable evolution of the (weakly perturbed) vortex soliton displayed
in panel (a2).

power of the anisotropic vortex solitons on LI 𝛼2 with fixed parameters families is investigated by means of the numerical solution of the
linearized equations for small perturbations, and verified by systematic
𝐷1 = 𝐷2 = 1, 𝛼1 = 2, 𝑔 = 𝛾 = 1, (50) simulations of the perturbed evolution. In particular, we have obtained
is considered. One observes that the corresponding solitons are stable families of stable elliptically shaped fundamental solitons and vortex
only for 𝛼2 ≥ 1.94, i.e., very close to the limit case corresponding to the ones, with 𝑆 = 1, in the 2D model with different LIs acting the in 𝑥 and
normal (non-fractional) diffraction. The intensity and phase patterns 𝑦 directions. Stable vortex solitons with topological charges 𝑆 = 1, 2, 3
for 𝛼2 = 1.94 are displayed in Figs. 11(a2, a3), respectively. It is are produced in the isotropic 2D model.
seen that the intensity is shaped as an oval “crater’’, and the phase To extend the analysis reported in this paper, a challenging direction
also exhibits the anisotropy. The corresponding dynamical behavior is is to implement the dual-LI scheme in the 3D geometry.
displayed in Figs. 11 (a4, a5, a6). Furthermore, the comparison of the
shapes of the stable anisotropic vortex in Figs. 11(a2, a4, a5), and the CRediT authorship contribution statement
long evolution history, reported in panel 11(a6), clearly demonstrate
that, whilst the anisotropic vortex remains stable against splitting, it Ming Zhong: Writing – original draft, Methodology, Investiga-
demonstrates steady rotation, i.e., an example of the robust coupling tion, Formal analysis, Conceptualization. Yong Chen: Writing – review
between the spin and orbital angular momenta, see a review [68]. Note & editing, Methodology, Formal analysis, Conceptualization. Zhenya
that the fundamental solitons in the same anisotropic model do not Yan: Writing – review & editing, Supervision, Methodology, Investiga-
exhibit rotation, as seen in Fig. 5(c1). tion, Funding acquisition, Formal analysis, Conceptualization. Boris A.
Malomed: Writing – review & editing, Methodology, Funding acquisi-
5. Conclusions and discussions tion, Formal analysis, Conceptualization.

In this work, we have proposed novel 1D and 2D models in the Declaration of competing interest
form of the dual-LI (dual-Lévy-index) FNLSEs (fractional nonlinear
Schrödinger equations) for the light propagation in fractional optical The authors declare that they have no known competing finan-
media. The models include two fractional-diffraction terms with dif- cial interests or personal relationships that could have appeared to
ferent values of the respective LIs. The nonlinearity, which provides influence the work reported in this paper.
for the stability of the model against the collapse, is taken as the
usual combination of self-focusing and defocusing cubic and quintic Data availability
terms, which is provided by various optical materials. The 2D set-
tings are elaborated in isotropic and anisotropic forms, the latter one
No data was used for the research described in the article.
corresponding to the diffraction operators with different LIs, acting
along two transverse coordinates. The models can be implemented in
optical waveguides with the effective fractional diffraction emulated by Acknowledgments
means of the properly patterned phase mask(s) acting on the Fourier-
transformed light fields [8,10]. Families of 1D and 2D fundamental This work was supported by the National Natural Science Foun-
and vortex solitons are obtained in the numerical form and by means dation of China (Grant No. 11925108) and Israel Science Foundation
of the VA (variational approximation). The stability of the soliton (Grant No. 1695/22).

12
M. Zhong et al. Physica D: Nonlinear Phenomena 462 (2024) 134157

References [35] S.N. Vlasov, V.A. Petrishchev, V.I. Talanov, Averaged description of wave beams
in linear and nonlinear media (the method of moments), Radiophys. Quantum
[1] N. Laskin, Fractional quantum mechanics, Phys. Rev. E 62 (2000) 3135. El. 14 (1971) 1062.
[2] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: A [36] M. Quiroga-Teixeiro, H. Michinel, Stable azimuthal stationary state in quintic
fractional dynamics approach, Phys. Rep. 339 (2000) 1–77. nonlinear optical media, J. Opt. Soc. Amer. B 14 (1997) 2004.
[37] B.A. Malomed, L.C. Crasovan, D. Mihalache, Stability of vortex solitons in the
[3] G.M. Zaslavsky, Chaos, fractional kinetics, and anomalous transport, Phys. Rep.
cubic-quintic model, Physica D 161 (2002) 187.
371 (2002) 461–580.
[38] D. Mihalache, D. Mazilu, L.C. Crasovan, B.A. Malomed, F. Lederer, Three-
[4] N. Laskin, Fractional Schrödinger equation, Phys. Rev. E 66 (2002) 056108.
dimensional spinning solitons in the cubic-quintic nonlinear medium, Phys. Rev.
[5] N. Laskin, Fractional Quantum Mechanics, World Scientific, Singapore, 2018.
E 61 (2000) 7142.
[6] B. Mandelbrot, The Fractal Geometry of Nature, W. H. Freeman, New York, 1982.
[39] A. Desyatnikov, A. Maimistov, B.A. Malomed, Three-dimensional spinning soli-
[7] M. Riesz, L’intégrale de Riemann–Liouville et le problème de Cauchy, Acta Math. tons in dispersive media with the cubic-quintic nonlinearity, Phys. Rev. E 61
81 (1949) 1. (2000) 3107.
[8] S. Longhi, Fractional Schrödinger equation in optics, Opt. Lett. 40 (2015) 1117. [40] D. Mihalache, D. Mazilu, L.C. Crasovan, I. Towers, A.V. Buryak, B.A. Malomed,
[9] H. Kasprzak, Differentiation of a noninteger order and its optical implementation, L. Torner, J.P. Torres, F. Lederer, Stable spinning optical solitons in three
Appl. Opt. 21 (1982) 3287. dimensions, Phys. Rev. Lett. 88 (2002) 073902.
[10] S. Liu, Y. Zhang, B.A. Malomed, E. Karimi, Experimental realisations of the [41] Z. Yan, V.V. Konotop, Exact solutions to three-dimensional generalized nonlinear
fractional Schrödinger equation in the temporal domain, Nature Commun. 14 Schrödinger equations with varying potential and nonlinearities, Phys. Rev. E 80
(2023) 222. (2009) 036607.
[11] B.A. Stickler, Potential condensed-matter realization of space-fractional quantum [42] B.A. Malomed, Multidimensional Solitons, AIP Publishing, Melville, NY, 2022.
mechanics: The one-dimensional Lévy crystal, Phys. Rev. E 88 (2013) 012120. [43] E.L. Falcão-Filho, C.B. de Araújo, G. Boudebs, H. Leblond, V. Skarka, Robust
[12] F. Pinsker, W. Bao, Y. Zhang, H. Ohadi, A. Dreismann, J. Baumberg, Fractional two-dimensional spatial solitons in liquid carbon disulfide, Phys. Rev. Lett. 110
quantum mechanics in polariton condensates with velocity-dependent mass, (2013) 013901.
Phys. Rev. B 92 (2015) 195310. [44] Reyna. A.S., C.B. de Araújo, High-order optical nonlinearities in plasmonic
[13] B.A. Malomed, Optical solitons and vortices in fractional media: A mini-review nanocomposites – A review, Adv. Opt. Phot. 9 (2017) 720–774.
of recent results, Photonics 8 (2021) 353. [45] S. Roy, F. Biancalana, Formation of quartic solitons and a localized continuum
[14] B.A. Malomed, Basic fractional nonlinear-wave models and solitons, Chaos 34 in silicon-based slot waveguides, Phys. Rev. A 87 (2013) 025801.
(2024) 022102. [46] K.K.K. Tam, T.J. Alexander, A. Blanco-Redondo, C.M. de Sterke, Generalized
[15] W.P. Zhong, M.R. Belić, B.A. Malomed, Y. Zhang, T. Huang, Spatiotemporal dispersion Kerr solitons, Phys. Rev. A 101 (2020) 043822.
accessible solitons in fractional dimensions, Phys. Rev. E 94 (2016) 012216. [47] A. Blanco-Redondo, C.M. De Sterke, J.E. Sipe, T.F. Krauss, B.J. Eggleton, C.
[16] C. Huang, L. Dong, Gap solitons in the nonlinear fractional Schrödinger equation Husko, Pure-quartic solitons, Nature Commun. 7 (2016) 10427.
with an optical lattice, Opt. Lett. 41 (2016) 5636. [48] A. Aceves, A. Copeland, Spatio-temporal dynamics in the mixed fractional
nonlinear Schrödinger equation, Front. Photon. 3 (2022) 977343.
[17] L. Zeng, J. Zeng, One-dimensional gap solitons in quintic and cubic-quintic
[49] F. Kh. Abdullaev, M. Salerno, Gap-Townes solitons and localized excitations in
fractional nonlinear Schrödinger equations with a periodically modulated linear
low-dimensional Bose–Einstein condensates in optical lattices, Phys. Rev. A 72
potential, Nonlinear Dynam. 98 (2019) 985.
(2005) 033617.
[18] M. Zhong, Z. Yan, Formation of multi-peak gap solitons and stable excitations
[50] L. Zeng, Y. Zhu, B.A. Malomed, D. Mihalache, Q. Wang, H. Long, Y. Cai, X.W. Lu,
for double-Lévy-index and mixed fractional NLS equations with optical lattice
J. Li, Quadratic fractional solitons, Chaos Solitons Fractals 154 (2022) 111586.
potentials, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 479 (2023) 20230222.
[51] Y. Qiu, B.A. Malomed, D. Mihalache, X. Zhu, L. Zhang, Y. He, Soliton dynamics
[19] L. Zeng, J. Zeng, One-dimensional solitons in fractional Schrödinger equation
in a fractional complex Ginzburg–Landau model, Chaos Solitons Fractals 131
with a spatially periodical modulated nonlinearity: nonlinear lattice, Opt. Lett.
(2020) 109471.
44 (2019) 2661. [52] M. Zhong, Two-dimensional fractional  -symmetric cubic-quintic NLS equa-
[20] Y. Qiu, B.A. Malomed, D. Mihalache, X. Zhu, X. Peng, Y. He, Stabilization of tion: Double-loop symmetry breaking bifurcations, ghost states and dynamics,
single and multi-peak solitons in the fractional nonlinear Schrödinger equation Physica D 448 (2023) 133727.
with a trapping potential, Chaos Solitons Fractals 140 (2020) 110222. [53] C. Klein, C. Sparber, P. Markowich, Numerical study of fractional nonlinear
[21] X. Yao, X. Liu, Off-site and on-site vortex solitons in space-fractional photonic Schrödinger equations, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 470
lattices, Opt. Lett. 43 (2018) 5749. (2014) 20140364.
[22] P. Li, B.A. Malomed, D. Mihalache, Vortex solitons in fractional nonlinear [54] J. Yang, Nonlinear Waves in Integrable and Nonintegrable Systems, SIAM, 2010.
Schrödinger equation with the cubic-quintic nonlinearity, Chaos Solitons Fractals [55] N.G. Vakhitov, A.A. Kolokolov, Stationary solutions of the wave equation in the
137 (2020) 109783. medium with nonlinearity saturation, Radiophys. Quantum Elect. 16 (1973) 783.
[23] Q. Wang, G. Liang, Vortex and cluster solitons in nonlocal nonlinear fractional [56] R.L. Pego, H.A. Warchall, Spectrally stable encapsulated vortices for nonlinear
Schrödinger equation, J. Opt. 22 (2020) 055501. Schrödinger equations, J. Nonlinear Sci. 12 (2002) 347.
[24] K. Li, P. Kevrekidis, D.J. Frantzeskakis, C.E.R. üter, D. Kip, Revisiting the- [57] D. Mihalache, D. Mazilu, I. Towers, B.A. Malomed, F. Lederer, Stable spatiotem-
symmetric trimer: Bifurcations, ghost states and associated dynamics, J. Phys. poral spinning solitons in a bimodal cubic-quintic medium, Phys. Rev. E 67
A: Math. Theor. 46 (2013) 375304. (2003) 056608.
[25] H. Susanto, R. Kusdiantara, N. Li, O.B. Kirikchi, D. Adzkiya, E.R.M. Putri, T. [58] L. Dong, F. Ye, J. Wang, T. Cai, Y.P. Li, Internal modes of localized optical vortex
Asfihani, Snakes and ghosts in a parity-time-symmetric chain of dimers, Phys. soliton in a cubic-quintic nonlinear medium, Physica D 194 (2004) 219.
Rev. E 97 (2018) 062204. [59] N. Dror, B.A. Malomed, Symmetric and asymmetric solitons and vortices in
[26] P. Li, B.A. Malomed, D. Mihalache, Symmetry-breaking bifurcations and ghost linearly coupled two-dimensional waveguides with the cubic-quintic nonlinearity,
states in the fractional nonlinear Schrödinger equation with a  -symmetric Physica D 240 (2011) 526.
potential, Opt. Lett. 46 (2021) 3267. [60] P. Li, B.A. Malomed, D. Mihalache, Metastable soliton necklaces supported
[27] P. Li, R. Li, C. Dai, Existence, symmetry breaking bifurcation and stability of by fractional diffraction and competing nonlinearities, Opt. Express 28 (2020)
two-dimensional optical solitons supported by fractional diffraction, Opt. Express 34472.
29 (2021) 3193. [61] K.I. Pushkarov, D.I. Pushkarov, I.V. Tomov, Self-action of light beams in
nonlinear media: soliton solutions, Opt. Quantum. Electron. 11 (1979) 471.
[28] M. Zhong, L. Wang, P. Li, Z. Yan, Spontaneous symmetry breaking and ghost
[62] Z.E. Birnbaum, B.A. Malomed, Families of spatial solitons in a two-channel
states supported by the fractional nonlinear Schr ödinger equation with focusing
waveguide with the cubic-quintic nonlinearity, Physica D 237 (2008) 3252.
saturable nonlinearity and  -symmetric potential, Chaos 33 (2023) 013106.
[63] R.L. Frank, E. Lenzmann, Uniqueness and nondegeneracy of ground states for
[29] M. Zhong, Z. Yan, Spontaneous symmetry breaking and ghost states in two-
(−𝛥)𝑠 𝑄 + 𝑄 − 𝑄𝛼+1 = 0 in R, 2010, arXiv:1009.4042.
dimensional fractional nonlinear media with non-Hermitian potential, Commun.
[64] M. Chen, S. Zeng, D. Lu, W. Hu, Q. Guo, Optical solitons, self-focusing, and wave
Phys. 6 (2023) 92.
collapse in a space-fractional Schrödinger equation with a Kerr-type nonlinearity,
[30] P. Li, H. Sakaguchi, L. Zeng, X. Zhu, D. Mihalache, B.A. Malomed, Second-
Phys. Rev. E 98 (2018) 022211.
harmonic generation in the system with fractional diffraction, Chaos Solitons [65] H. Sakaguchi, B. Li, B.A. Malomed, Creation of two-dimensional composite
Fractals 173 (2023) 113701. solitons in spin–orbit-coupled self-attractive Bose–Einstein condensates in free
[31] S. Kumar, P. Li, B.A. Malomed, Domain walls in fractional media, Phys. Rev. E space, Phys. Rev. E 89 (2014) 032920.
106 (2022) 054207. [66] J. Qin, G. Dong, B.A. Malomed, Stable giant vortex annuli in microwave-coupled
[32] L. Bergé, Wave collapse in physics: Principles and applications to light and atomic condensates, Phys. Rev. A 94 (2016) 053611.
plasma waves, Phys. Rep. 303 (1998) 259. [67] B.A. Malomed, Two-dimensional solitons in nonlocal media: A brief review,
[33] G. Fibich, The Nonlinear Schrödinger Equation: Singular Solutions and Optical Symmetry 14 (2022) 1565.
Collapse, Springer, New York, 2015. [68] K.Y. Bliokh, F.J. Rodríguez-Fortuño, F. Nori, A.V. Zayats, Spin–orbit interactions
[34] E. Garmire, R.Y. Chiao, C.H. Townes, Dynamics and characteristics of the of light, Nature Photon. 9 (2015) 796–808.
self-trapping of intense light beams, Phys. Rev. Lett. 16 (1966) 347.

13