Introduction

Solitons, alias localized waves, are stable objects balanced by diffraction/dispersion and nonlinearity1,2,3. The stable propagation and self-trapping of two- and three-dimensional (2D and 3D) solitons have recently received reinvigorating scientific interests in physics and beyond4. A fundamental challenge in this field of research is the stabilization of solitons in multidimensional coordinates, since the 2D and 3D solitons in free space are always unstable and undergo respectively critical and supercritical collapses arising from catastrophic self-focusing nonlinearity5,6. It is therefore intriguing to suppress these collapses in the recent soliton research7,8.

The stabilization of multidimensional localized states usually relies on nonuniform media. A common way is to introduce linear periodic potentials (alias linear lattices)7,8,9,10 to a material with uniform nonlinearity, e.g., a periodic transverse modulation of the refractive index (photonic crystal and lattice) in optics11,12 or an optical lattice generated by counter-propagating laser beams in Bose–Einstein condensates13,14,15,16,17. The periodic potentials play a paramount stabilizing part for D-dimensional localized modes including fundamental solitons and more complicated states like gap solitons, multihump states, and vortex solitons4,7,8. Another potential scheme is assisted by nonlinear pseudopotentials—the well-known nonlinear lattices7,18,19,20. Although nonlinear lattices with smooth variation of nonlinearity may support various species of solitons, to stabilize multidimensional solitons using sinusoidal nonlinear lattices remains an open area of research7. Further, the combined linear-nonlinear lattices yield stabilization of various soliton families21,22,23.

Conversely, materials which exhibit different orders of nonlinearity, including the well-known cubic–quintic model with competing nonlinearities (self-focusing cubic and self-defocusing quintic nonlinear terms), can create stable multidimensional solitons and suppress the above-mentioned critical- and supercritical collapses4,7,8,21,22. Experimentally, fundamental solitons and vortex solitons in 2D and 3D geometries have been observed in CS\({}_{2}\)24,25 and in suspensions of metallic nanoparticles26. Theoretical predictions27,28 and experimental observations29,30,31,32 very recently confirmed the existence of stable multidimensional states of ultracold bosonic gases8, appearing as the so-called quantum droplets27,28, where the Bose–Bose mixture collapses are suppressed by quantum fluctuations33.

Beyond strategies affording the stabilization of multidimensional solitons in standard quantum mechanics, two extensions of which have been elaborated in past decades. One is the non-Hermitian system whose Hamiltonian is parity-time (\({\mathcal{P}}{\mathcal{T}}\)) symmetric and, due to this unique property, such system—counterintuitively—exhibits entirely real eigenvalue spectra34,35,36. In recent years, the studies of various types of solitons37 and their propagation dynamics governed by different \({\mathcal{P}}{\mathcal{T}}\)-symmetric systems integrated with \({\mathcal{P}}{\mathcal{T}}\) lattices (periodic potentials with alternating gain and loss regions) are fruitful38,39,40,41,42. Another extension, made by Nick Laskin, is space-fractional quantum mechanics (SFQM)43,44. The SFQM is achieved if the Lévy flights replaced the Brownian trajectories in Feynman path integrals, leading to a fractional diffraction order characterized by Lévy index \(\alpha\) (\(1\, < \, \alpha \leq 2\)) rather than the standard one with \(\alpha =2\). As the Schrödinger equation roots in the path integral over Brownian paths, the path integral over Lévy trajectories results into fractional Schrödinger equation (FSE)43,44,45.

FSE is currently attracting a great deal of attention in different branches of physics, including condensed-matter physics46,47,48 and quantum physics49,50,51. Importantly, a new hallmark of FSE development published in 2015 by Longhi52 found that a fractional quantum harmonic oscillator can be realized based on transverse laser beam dynamics in aspherical optical cavities, pushing the studies of fractional physical system into optics. Such work triggers a growing prosperity of fractional models for digging deep into the propagation properties of the possible beam solutions in both the linear and nonlinear systems, striking examples include: Gaussian beams propagation53,54, conical diffraction of light (\({\mathcal{PT}}\) symmetry) in FSE with a periodic \({\mathcal{PT}}\)-symmetric potential55, accessible solitons in FSE with a harmonic potential56, solitons57,58, and modulational instability59 in purely nonlinear fractional Schrödinger equation (NLFSE), and diverse soliton families in the NLFSE setting by adding to its linear60,61,62,63 or nonlinear20 part with a periodic potential (the aforementioned optical and nonlinear lattices, respectively). Although the NLFSE provides a fertile ground for exploring various types of solitons, their existence and stability property supported by the cubic–quintic nonlinearities and 2D linear periodic potential (optical lattice) are yet to be revealed.

In this work, we propose and demonstrate, theoretically and numerically, a framework of 2D NLFSE—describing light propagation in a nonlinear periodic system with an optical lattice and attractive (self-focusing)-repulsive (self-defocusing) cubic–quintic nonlinearities—which can suppress the critical collapse mentioned before. We reveal that the model produces a variety of stable soliton families, including 2D fundamental gap and vortical solitons as well as gap soliton clusters (solitons are always unstable in the quintic-only model). Noticeably, the Lévy index \(\alpha\) and propagation constant \(b\) impact deeply upon the profiles of soliton families, the side peaks of which with higher modulations when decreasing \(\alpha\) or increasing \(b\). By performing the linear-stability analysis and direct numerical simulations, the stability and instability diagrams of all the soliton families are acquired. A detailed insight into the dynamic properties of solitons further shows that the solitons are robustly stable in the middle of the band gaps of the underlying linear Bloch spectrum, while unstable near the edges of the band gaps; and the stability of the solitons is moderately influenced by nonlinear strength. Our work proposes a feasible scheme to stabilize 2D localized modes against critical collapse by considering the fractional diffraction order to light propagation in periodic physical systems with competing self-focusing and self-defocusing nonlinearities in cubic–quintic nonlinear terms, thus offering a new avenue to investigate the existence and dynamic properties of 2D localized modes by managing the diffraction order and tunable band gaps of the systems.

Results

Theoretical model and its linear spectrum

We start from the cubic–quintic NLFSE (with competing focusing-defocusing nonlinearities), which describes a laser beam propagation in a nonlinear medium and can be expressed in scaled form:

$$i\frac{\partial E}{\partial z}=\frac{1}{2}{\nabla }^{\alpha }E+V(x,y)E-\gamma {\left|E\right|}^{2}E+{\left|E\right|}^{4}E,$$
(1)

where \(E\) and \(z\) represent the field amplitude and propagation distance, respectively, \(\gamma \, > \, 0\) being the coefficient of self-focusing cubic nonlinearity. Note that fractional Laplacian \({\nabla }^{\alpha }={(-{\partial }_{x}^{2})}^{\alpha /2}+{(-{\partial }_{y}^{2})}^{\alpha /2}\) acts on two transverse coordinates \(x\) and \(y\), with \(\alpha\) (\(1 \, < \, \alpha \leq 2\)) corresponds to Lévy index. In particular, Eq. (1) restores to the conventional nonlinear Schrödinger equation at \(\alpha =2\). The remaining parameter \(V\) is the linear (square) potential trap, which yields

$$V={V}_{0}\left[{{\rm{sin}}}^{2}(x)+{{\rm{sin}}}^{2}(y)\right],$$
(2)

where \({V}_{0}\, > \, 0\) stands for the amplitude (modulation depth) of optical lattice.

To explore the existence of 2D nonlinear localized gap modes, we should first give out the band gap structure of guided Bloch modes (linear Bloch waves) stemming from the linearized model of Eq. (1) (see “Methods”). For the given expression of the sinusoidal square optical lattice (2), whose contour plot is portrayed in Fig. 1a, the corresponding band gap structure of Bloch-wave spectrum, as mentioned above, can be described and found by adopting the analogies with the solid energy band theory of crystalline materials and which, in the way as usual, should be based on the lowest (first) Brillouin zone of the lattice potentials. Such reduced Brillouin zone of the square lattice under considered is displayed in Fig. 1b. In the view of these foundations, we have plotted the dispersion diagram for guided Bloch waves with the variation of Lévy index \(\alpha\) in Fig. 1c, and shown such dispersion spectrum for a particular value of \(\alpha\) (\(\alpha =1.3\)) in Fig. 1d. It is observed from the former panel that both the first and second finite band gaps widen with an increasing \(\alpha\), conforming to the previous reported results in 1D situations63. As for the latter panel, choosing \(\alpha =1.3\) is qualified to guarantee the moderate wideness of the first two band gaps, within which various localized gap modes may be stabilized by the competing cubic–quintic nonlinear terms.

Fig. 1: Shape, Brillouin zone, and linear spectrum of the 2D linear (index) potential.
figure 1

a Contour plot of a sinusoidal square lattice in Cartesian space, blue shading corresponds to potential minima. b The first reduced Brillouin zone of such 2D square lattice in reciprocal lattice space; tagged are the high-symmetry points. c Band gap structure for guided Bloch waves with different \(\alpha\). d The dispersion spectrum at Lévy index \(\alpha =1.3\), with the colored curves denoting Bloch bands. \({V}_{0}=6\) and \(\gamma =0.5\) are used throughout this paper.

Full size image

Formation of two-dimensional nonlinear localized gap modes

Having obtained the finite band gaps of the underlying linear spectrum to the model (Eq. (1)), in the next part of this research work we are going to explore the formation of 2D nonlinear localized gap modes (with their corresponding propagation constants lying in such band gaps) supported by the interplay of the given linear lattice and the competing cubic–quintic nonlinearities. Several classes of soliton families are identified, such as fundamental gap solitons and vortical ones as well as the gap-type soliton clusters. The physical mechanism and theoretical analysis are also given out to all the localized modes thus found.

To understand how the Lévy index (\(\alpha\)) in the fractional Laplacian (diffraction) term affects the shape of fundamental gap solitons, we plotted the spatial structure of such fundamental localized modes at two different values of \(\alpha\) by means of 3D plots, top views and contour profiles in Fig. 2a–f. Comparing the solitons’ shape in these panels, we can see that the fundamental gap soliton is in a wavy shape at smaller \(\alpha\) (Fig. 2a–c), while such wavy modulation of the gap soliton disappears at relatively large \(\alpha\) (Fig. 2d–f); this characteristic resembles those found in analog soliton-based systems60,63.

Fig. 2: Typical examples of fundamental gap solitons at different values of Lévy index.
figure 2

a–c \(\alpha =1.1\), d–f \(\alpha =1.9\), propagation constant b = −3.3 for all. The 3D profiles (a, d), top views (b, e), and contour profiles (c, f) are shown for the localized gap modes.

Full size image

Our systematic simulations by exploiting the undermentioned linear stability and direct perturbed propagation methods (see “Methods” for detailed descriptions) have demonstrated the presence of such fundamental gap solitons dwelling in the first and second band gaps, and identified their stability and instability properties, which are accumulated in the associated \(P\)\(\, \sim \,\)\(b\) curve in Fig. 3a. From such curve one can clearly observe that the stability region is much broader when preparing the gap solitons in the first band gap than that in the second gap; and intuitively, the unstable gap solitons emerge once the propagation constant is near the edge of the finite band gap, as within where the effective mass of the gap solitons may change the sign7. Note that such curves (\(P\)\(\, \sim \,\)\(b\)) for multipole gap solitons (or gap soliton clusters) with numbers of peaks \(n=4\) and \(n=16\) are showed respectively in Fig. 3b, c.

Fig. 3: Soliton power \(P\) versus propagation constant \(b\) for localized gap modes.
figure 3

The dependence \(P(b)\) with different numbers of peaks \(n\): (a) fundamental gap soliton (\(n=1\)); multipole gap solitons with \(n=4\) (b), and \(n=16\) (c). The blue and red lines stand for stable and unstable domains. For the sake of discussion, profiles of solitons marked by circles (A1, A2, A3), (B1, B2, B3), and (C1, C2, C3) will be displayed in the following figures.

Full size image

Three representative scenarios of such fundamental gap solitons, residing on the first two band gaps, are severally depicted in Fig. 4(a, d), (b, e) and (c, f) (which correspond to the marked circles (A1, A2, A3) in Fig. 3a). When placing in much higher band gap, e.g., the second band gap, the gap soliton becomes more localized than that in the first band gap, as comparing the Fig. 4a, d with 4b, e. It is observed from the Fig. 4c, f that the gap soliton becomes highly confined mode occupying several lattice sites when the propagation constant \(b\) is close to the edge of the second band gap, with an emergence of a cusplike behavior, such feature is typical for the gap modes in other physical settings7,16. The eigenvalues, calculated from the linear stability equations (see the Eq. (4) in “Methods”), of these three gap modes are displayed in Fig. 4g–i, showing that the former two modes are linearly stable while unstable for the latter mode which lies inside the upper edge of the second band gap.

Fig. 4: Characteristic shapes and eigenvalues of fundamental gap solitons.
figure 4

Contour plots at different values of propagation constant b: (a) b = −3.5; (b) b = −5.1; (c) b = −5.7. Subplots (d–f) and (g–i) are the corresponding profiles and eigenvalues of fundamental gap solitons displayed in subplots (a–c).

Full size image

On the other hand, our research model (Eq. (1)) upholds families of higher-order solitons, which we call gap soliton clusters. Depicted in Fig. 5 are the quadruple-mode solitons of such soliton clusters, with a distance (\(\bigtriangleup s\)) between isolated gap solitons doubling to the period of the lattice. Because of the existence of repelling interaction between the each adjacent solitons of the soliton clusters, the stability diagram of the quadruple-mode solitons, portrayed as the dependency \(P(b)\) in Fig. 3b, is much narrower than that for the fundamental gap solitons. The stability and instability properties of such multiple solitons with the number of isolated gap solitons \(n=4\) are further validated by the corresponding eigenvalues analysis, which are showed in the bottom of Fig. 5.

Fig. 5: Contour plots, profiles, and eigenvalues of multiple solitons with \(n=4\).
figure 5

Multiple gap solitons at different values of b: (a) b = −3.5; (b) b = −5.2; (c) b = −5.7. Subplots (d–f) and (g–i) are the corresponding profiles and eigenvalues of multiple gap solitons displayed in subplots (a–c).

Full size image

The successful stabilization of quadruple-mode solitons enlightens us to generate multiple solitons with more number of gap solitons, which can be naturally achieved by the model under studied. For instance, the gap soliton clusters with the neighboring distance \(\bigtriangleup s=2\pi\) and \(n=16\) solitons arranging as \(4\times 4\) configuration, are displayed as stable modes localized in the first and second band gaps in Fig. 6a, d, b, e, respectively, and as unstable object when its propagation constant approaches the edge of the band gap in Fig. 6c, f. We stress that the multiple gap solitons or gap soliton clusters reported here are equivalent to the recently found truncated nonlinear Bloch waves or gap waves in nonlinear optical media and ultracold atoms17,23,64,65,66,67,68.

Fig. 6: Contour plots, profiles, and eigenvalues of multiple solitons with \(n=16\).
figure 6

Multiple gap solitons at different values of b: (a) b = −3.5; (b) b = −5.4; (c) b = −5.7. Subplots (d–f) and (g–i) are the corresponding profiles and eigenvalues of these 16-mode clusters displayed in subplots (a–c).

Full size image

Stable propagations of the localized gap modes thus found (fundamental gap solitons and soliton clusters), located inside the first and second band gaps, are depicted respectively in Fig. 7 which, in practice, demonstrates the stability properties (typical for solitons) of the localized gap modes in a direct and vivid way. Typical propagation properties of unstable gap modes are depicted in Fig. 8, indicating that the unstable modes exhibit a weakly instability (whereas the main peaks can sustain, their side peaks grow quickly during the propagation), which is in an excellent agreement with the relevant linear eigenvalues results as shown in the bottom lines of the Figs. 4–6.

Fig. 7: Stable propagations of localized gap modes.
figure 7

Fundamental gap solitons at b = −3.5 (a) and b = −5.1 (d); multiple gap solitons with number of peaks \(n=4\) at b = −3.5 (b) and b = −5.2 (e); multiple gap solitons with \(n=16\) at b = −3.5 (c) and b = −5.4 (f).

Full size image
Fig. 8: Unstable propagations of localized gap modes at b = −5.7.
figure 8

Fundamental gap soliton (a); multiple gap solitons with \(n=4\) (b), and \(n=16\) (c).

Full size image

An interesting result is the possibility to find stable localized gap vortices, the gap solitons carrying with topological charge m. A typical gap vortex is usually constructed as quadruple-mode bound states of identical solitons with modulated phase structure. Characteristic examples of such gap vortices with m = 1 and spacing \(\bigtriangleup s=2\pi\) are shown in Fig. 9. It is shown that the complicate phase structures of these gap vortices are modulated spatially owning to the influence of the lattice potential. Eigenvalues of these vortical solutions, depicted in the bottom of Fig. 9, indicate that there is weak linear instability for the unstable gap vortices, such effect equates to that for the localized gap modes without vorticity (m = 0) (see the panels in the bottom right corner of Figs. 4–6) and is verified as well in their direct perturbed propagation in Fig. 10. For defined values of propagation b and cubic nonlinear coefficient \(\gamma\), we find that the stability and instability regions of the vortex gap solitons alter with the variation of Lévy index α , as attested by the dependence P(α) in Fig. 10a.

Fig. 9: Contour plots, phases, and eigenvalues of vortex gap solitons appearing as quadruple-mode bound states.
figure 9

Vortex gap solitons with vortex charge m = 1 and b = −4, and at different values of Lévy index: (a) α  =  1.2; (b) α  =  1.2; (c) α  = 1.9. The phase structures and eigenvalues of vortex modes in subplots (a–c) are, respectively, displayed in subplots (d–f) and (g–i).

Full size image
Fig. 10: Soliton power P versus α, and propagation dynamics of vortex gap solitons.
figure 10

a Dependence P (b) of vortex gap solitons (quadruple-mode vortices with n = 4) at vortex charge m = 1 and b = −4. The blue and red lines in a correspond to stable and unstable domains. Propagations of vortex gap solitons shown as marked circles (D1, D2, D3) in a: unstable gap vortexes at α = 1.2 (b) and α  =  1.3 (c); (d) stable gap vortex at α  =  1.9.

Full size image

Discussion

As previously described, the 2D solitons created in pure self-focusing settings (uniform Kerr media) undergo critical collapse, and their vortex counterparts suffer a greater azimuthal instability. To suppress such adverse factors, scientists usually resort to two alternative methods—the cubic–quintic model with competing self-focusing and self-defocusing nonlinearities, and the introduction of external linear potentials such as optical lattices. By combing such two methods, we here numerically explored and with an analysis of the generation of localized modes in the 2D cubic–quintic nonlinear FSE with an optical lattice. We discovered that such model gives rise to a class of 2D stable soliton families, which include 2D fundamental gap solitons and their clusters as well as vortex solitons. Like previous studies in similar but different models, we demonstrated that all the gap soliton families are always unstable physical objects in the quintic-only model (the cubic–quintic model discarding the cubic term). We emphasize that, the two physical parameters, the Lévy index α and propagation constant b, impose a nontrivial effect upon the profiles of soliton families of all kinds, the solitons accompanied by a higher modulated side peaks if we simply minish α or enlarge b. We identify the stability and instability diagrams of all the soliton families by right of the linear-stability analysis by calculating the eigenvalues of the underlying soliton solutions against linear infinitesimal background perturbations, and the direct numerical simulations. We conclude that the solitons are stable within the band gaps (of the underlying linear system) and unstable close to the band edges.

An obvious conceptional extension of our study is to reveal the underlying physical mechanism for critical collapse related to physical parameters including Lévy index α (diffraction order) and order of a single nonlinearity (e.g., cubic only) as well as the strength of optical lattices. Rich localized modes and nonlinear dynamics are yet to be explored in an incommensurate structure—combined linear and nonlinear lattices21,22,23. Another interesting issue is to study the localization and delocalization properties of waves in other periodic lattices69 with regolabile (fractional) diffraction.

Methods

Additional details on the numerical procedure

To obtain the stationary solution U, we first rewrite the field amplitudes as E = U exp(ibz) with the propagation constant \(b\) and plug them into the Eqs. (1), then we get the following stationary equation:

$$-bU=\frac{1}{2}{\nabla }^{\alpha }U+V(x,y)U-\gamma {\left|U\right|}^{2}U+{\left|U\right|}^{4}U.$$
(3)

An essential physical quantity is the soliton power P, which is defined as \(P=\int \int {\left|U(x,y)\right|}^{2}dxdy\).

For another, the linear stability analysis of the thus-found gap solitons is a critical issue. To this, we express the perturbed field amplitude as \(E=[U(x,y)+p(x,y){\mathrm{exp}}(\lambda z)+{q}^{* }(x,y){\mathrm{exp}}({\lambda }^{* }z)]{\mathrm{exp}}(ibz)\), here \(U(x,y)\) represents undisturbed field amplitude, p(x, y) and \({q}^{* }(x,y)\) being small perturbations at eigenvalue \(\lambda\). Substituting such perturbed solution into the Eq. (1), we readily obtain the following eigenvalue problem:

$$\left\{\begin{array}{lll}i\lambda p= \ +\,\frac{1}{2}{\nabla }^{\alpha }p+[b+V(x,y)]p \hfill \\ -\, \gamma U(2{U}^{* }p+Uq)+{U}^{2}{U}^{* }(3{U}^{* }p+2Uq),\hfill \\ i\lambda q= -\,\frac{1}{2}{\nabla }^{\alpha }q-[b+V(x,y)]q \hfill \\ +\,\gamma {U}^{* }(2Uq+{U}^{* }p)-{({U}^{* })}^{2}U(3Uq+2{U}^{* }p).\end{array}\right.$$
(4)

We stress that from the above eigenvalue Eq. (4), the perturbed solutions are stable as long as all the real parts \(({\lambda }_{{\rm{R}}})\) are zero (\({\lambda }_{{\rm{R}}}=0\)). The numerical recipes we used follow the ways: the calculations of stationary solutions of Eq. (3) and numerical propagation simulations of Eq. (1) under small initial perturbations were, respectively, performed based on the modified squared-operator method and split-step Fourier method70.